Probabilistic forecasts for the 2026 FIFA World Cup are obtained by using a hybrid model that combines data, expert insights, and advanced statistical models. The favorite is Spain, closely followed by England, France, and Germany.

Football fans around the world are looking forward to the kick-off of the 2026 FIFA World Cup in Canada, Mexico, and the United States next week. 48 of the best teams from all around the world will compete from 11 June to 19 July to determine the new World Champion. In anticipation of the tournament the big question is who among the teams will succeed, who will drop out, and who will eventually prevail. While it is, of course, not yet possible to give definitive answers to these questions, we are able to provide probabilistic forecasts for all possible matches using a refined machine learning algorithm. This allows us to explore the likely course of the tournament by simulation.

Winning probabilities

The forecast is based on a machine learning algorithm that blends a variety of different sources of information: An ability estimate for every team based on historic matches; an ability estimate for every team based on odds from 24 bookmakers; average ratings of the players in each team based on their individual performances in their home clubs and national teams; the average market value of all players in each team according to a wisdom-of-the-crowd approach; further team and country covariates (e.g., FIFA and Elo ratings or GDP). A machine learning algorithm is trained on the results of all major football tournaments (Men’s World Cups and Euros) between 2006 and 2024 and then applied to current information to obtain a forecast for the 2026 FIFA World Cup. More specifically, the algorithm estimates the predicted number of goals for all possible matches between all 48 teams in the tournament. Based on the predicted goals the probabilities for each potential outcome (i.e., 0-0, 1-0, 0-1, 2-0, etc.) in each of these matches can be computed from a bivariate Poisson distribution (here: assuming independence). This allows us to simulate all matches in the group phase and which teams proceed to the knockout stage and who eventually wins. Repeating the simulation 100,000 times yields winning probabilities for each team. The results show that Spain is the favorite for the title with a winning probability of 14.5%, closely followed by England and France, both with 12.4%, and Germany with 11.2%. The winning probabilities for all teams are shown in the barchart below with more information linked in the interactive full-width version.

Interactive full-width graphic

The study has been conducted by an international team of researchers: Andreas Groll, Agamyrat Hanekov, Lars Magnus Hvattum, Rouven Michels, Gunther Schauberger, Elina Sukhanova, Sebastian Witte, Achim Zeileis. The basic idea for the forecast is to proceed in two steps. In the first step, sophisticated statistical models as well as expert insights are employed to determine the strengths of all teams and their players using disparate sets of information. In the second step, a machine learning algorithm decides how to best combine the strength estimates with other information about the teams.

• Historic information: Match abilities.
  An ability estimate is obtained for every team based on “retrospective” data, namely all historic national matches over the last 8 years (freely curated by Mart Jürisoo on Kaggle). A bivariate Poisson model with team-specific fixed effects and assuming independence is fitted to the number of goals scored by both teams in each match. However, rather than equally weighting all matches to obtain average team abilities (or team strengths) over the entire history period, an exponential weighting scheme is employed. This assigns more weight to more recent results and thus yields an estimate of current team abilities. More details can be found in Ley, Van de Wiele, Van Eetvelde (2019).

• Future expectation: Bookmaker consensus abilities.
  Another ability estimate for every team is obtained based on “prospective” data, namely the odds of 24 international bookmakers that reflect their expert expectations for the tournament. Using the bookmaker consensus model of Leitner, Zeileis, Hornik (2010), the bookmaker odds are first adjusted for the bookmakers’ profit margins (“overround”) and then averaged (on a logit scale) to obtain a consensus for the winning probability of each team. To adjust for the effects of the tournament draw (that might have led to easier or harder groups for some teams), an “inverse” simulation approach is used to infer which team abilities are most likely to lead up to the consensus winning probabilities.

• Individual player contributions: Average player ratings.
  To infer the “contributions of individual players” in a match, the plus-minus player ratings of Pantuso & Hvattum (2021) dissect all matches with a certain player (both on club and on national level) into segments, e.g., between substitutions. Subsequently, the goal difference achieved in these segments is linked to the presence of the individual players during that segment. This yields individual ratings for all players that can be aggregated to average player ratings for each team.

• Wisdom of the crowd: Average market values:
  Another way to reflect the current quality and the future potential of each player in a team is to consider their expected market value. As the real market values are unknown, the Transfermarkt web portal employs a “wisdom-of-the-crowd” approach to determine current expected market values for all players. These are based on discussions relying on publicly available data among the online community members of the portal and moderated and consolidated by expert community members and the portal’s employees.

• Combination with present status: Hybrid random forests.
  Finally, machine learning is used to combine these four highly aggregated and informative variables with a broad range of further relevant covariates reflecting the current states of the different teams and the countries they come from. Such a hybrid approach was first suggested by Groll, Ley, Schauberger, Van Eetvelde (2019). A random forest algorithm is trained to decide how to blend the different ability estimates with team-specific features that are typically less informative but still powerful enough to enhance the forecasts. The features considered comprise team-specific details (e.g., FIFA rank, Elo rating, number of Champions League players) as well as country-specifc socio-economic factors (such as GDP per capita). By combining a large ensemble of rather weakly informative regression trees in a random forest, the relative importances of all the covariates can be inferred automatically. The resulting predicted number of goals for each team can then finally be used to simulate the entire tournament 100,000 times.

Match probabilities

Using the forecasts from the machine learning algorithm yields the predicted number of goals for both teams in each possible match. The explanatory information used for this is the difference between the two teams in each of the variables listed above, i.e., the difference in historic match abilities (on a log scale), the difference in bookmaker consensus abilities (on a log scale), difference in average player ratings of the teams, difference in log market values, etc. The predicted number of goals for the two teams in each match can then be plugged as expectations into two independent Poisson distributions, from which we can compute the probability that a certain match ends in a win, a draw, or a loss. The same can be repeated in overtime, if necessary, and a coin flip is used to decide penalties, if needed.

The following heatmap shows for each possible combination of teams the probability that one team beats the other team in a knockout match. The color scheme uses green vs. purple to signal probabilities above vs. below 50%, respectively. The tooltips for each match in the interactive version of the graphic also print the probabilities for the match to end in a win, draw, or loss after normal time.

Interactive full-width graphic

Performance throughout the tournament

As the goals for both teams in every single match can be simulated with the approach described above, it is also straightfoward to simulate the entire tournament (here: 100,000 times) providing “survival” probabilities for each team across the different stages.

Interactive full-width graphic

Odds and ends

All our forecasts are probabilistic, clearly below 100%, and hence by no means certain. Although we can quantify this uncertainty in terms of probabilities from a multiverse of potential tournaments, it is far from being predetermined which of these potential tournaments we will eventually see during the actual tournament.

Nevertheless the probabilistic view provides us with some interesting insights: For example, compared to predictions for previous tournaments (see e.g., 2018, 2022), it is even more uncertain who will win the title as there are a number of teams with good (albeit none with very high) chances of winning the tournament. An important factor for this is the substantially increased size of the tournament with 48 teams (rather than the previous 32) and an additional knockout round. Also, the tournament draw is much more variable, because 8 of the 12 third-ranked teams proceed to the knockout stage with 495 (!) possible permutations for mapping groups to matches in the round of 32.

Moreover, comparing our forecasts to those based only on the bookmakers odds, it is striking that Germany is ranked 4th, closely behind the three top teams, while it is only ranked 7th by many bookmakers. Conversely, Brazil and Argentina are typically ranked higher by the bookmakers but perform worse in our machine-learning-calibrated simulation.

In any case, all of this means that the probabilistic forecasts leave a lot of room for surprises and excitement during the 2026 FIFA World Cup. But what is absolutely certain is that we look forward to an entertaining tournament as football fans (much more than as professional forecasters).

1. Introduction & Theoretical Framework

In modern electronic trading markets, algorithmic execution engines drive the vast majority of institutional order flows. Evaluating whether these independent, learning-driven trading algorithms behave competitively or tacitly coordinate has become a critical challenge for quantitative compliance, market microstructure design, and risk management.

This technical article implements an automated Strategic Audit Engine designed to evaluate algorithmic execution regimes in the Silver futures market (SI=F). Our framework is explicitly built upon the empirical and theoretical foundations laid out by Koulouris & Campajola (2026) in their groundbreaking paper, “Memory-Induced Supra-Competitive Outcomes Between Deep Reinforcement Learning Agents in Optimal Trade Execution” (arXiv:2605.20348v1, May 2026).

The Core Thesis: Supra-Competitive Outcomes via Memory Paths

Traditional regulatory frameworks look for explicit collusion (active communication or cartel setups). However, Koulouris & Campajola demonstrate a far more subtle phenomenon: when independent Deep Reinforcement Learning (DRL) agents are equipped with memory—meaning they learn from rolling windows of historical price trajectories—they naturally converge toward supra-competitive outcomes. These are states where joint rewards remain artificially high, or execution parameters naturally align to mimic cooperation, without any explicit information exchange.

To audit this behavior empirically, our engine models a symmetric duopoly market interaction. It maps the actual market execution path against two fundamental game-theoretic baselines:

• The Cooperative Boundary (TWAP / Pareto Frontier): An idealized, optimal trade execution path where volume is distributed evenly across time to minimize joint market impact and maximize long-term mutual utility.
• The Competitive Boundary (Nash Equilibrium): The aggressive, non-cooperative state where individual agents structurally undercut each other, driving execution shortfall parameters to their maximum baseline.

2. Technical Stack & Environmental Setup

To build a production-grade, reproducible multi-agent simulation pipeline, we leverage a hybrid data-science and deep-learning toolkit within the R ecosystem:

• tidyquant & tidyverse: Serve as our core data engineering layer, managing financial API queries, formatting continuous return matrices, and handling functional list columns.
• keras & tensorflow: Form the algorithmic backbone, allowing us to build, train, and run simultaneous forward/backward passes on Deep Q-Networks.
• ggtext & glue: Empower our visualization suite to parse inline HTML canvas rendering and handle dynamic string interpolations smoothly.

# 1. ENVIRONMENT SETUP
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyquant, tidyverse, ggtext, glue, keras, tensorflow)

3. Building the Double Deep Q-Network Topology

Following the paper’s thesis on symmetric duopoly interactions, we construct two structurally identical execution agents: agent_A and agent_B. Both utilize a Dense Neural Network (Multilayer Perceptron) architecture to approximate the action-value space, denoted as Q(s, a).

The state space contains 3 features: Price Deviation, Asset Volatility (sigma), and Relative Time Horizon. The output layer projects to 3 discrete strategic action coordinates via a linear activation function.

# 2. SYMMETRIC AGENT ARCHITECTURE
build_strategic_agent <- function(state_size = 3, action_size = 3) {
  model <- keras_model_sequential() %>%
    layer_dense(units = 32, activation = "relu", input_shape = c(state_size)) %>%
    layer_dense(units = 32, activation = "relu") %>%
    layer_dense(units = action_size, activation = "linear")
  
  model %>% compile(
    optimizer = optimizer_adam(learning_rate = 0.001),
    loss = "mse"
  )
  return(model)
}

# Initialize the competing agents
agent_A <- build_strategic_agent()
agent_B <- build_strategic_agent()

4. Parameterization & Historical Replay Buffer Ingestion

To anchor our agents in empirical reality, we pull 2 years of continuous daily settlement prices for Silver futures (SI=F). We define our microstructural bounds—such as the risk aversion parameter (gamma) and the permanent market impact vector (eta)—alongside a fixed strategic execution memory window (T = 10).

# 3. STRATEGIC PARAMETERS
T_horizon <- 10      # Strategic episode length (Memory window)
gamma_param <- 0.0001 # Risk aversion
eta_param <- 0.0005   # Market impact

# 4. HISTORICAL REPLAY DATA (2-Year Training Set)
silver_full <- tq_get("SI=F", from = Sys.Date() - 730) %>%
  filter(!is.na(close)) %>%
  mutate(returns = close / lag(close) - 1) %>%
  drop_na()

# Recent window for the final audit visualization
silver_recent <- tail(silver_full, T_horizon)

5. Dynamic Volatility Corridors

Rather than mapping market behavior against static thresholds, the audit engine computes a volatility-adaptive safety corridor. The boundaries dynamically expand and contract based on the asset’s realized standard deviation (sigma), isolating pure structural noise from intentional strategic maneuvers.

# 5. DYNAMIC SIGMA CORRIDORS
current_sigma <- sd(silver_recent$returns, na.rm = TRUE)
if(is.na(current_sigma)) current_sigma <- 0.01 

analysis_data <- silver_recent %>%
  mutate(
    twap_slope = current_sigma * 1.5, 
    nash_slope = current_sigma * 4.0,
    twap_path = first(close) * (1 - seq(0, first(twap_slope), length.out = n())),
    nash_path = first(close) * (1 - seq(0, first(nash_slope), length.out = n())),
    lower_safety_limit = nash_path * (1 - current_sigma)
  )

6. The Joint Training Replay Engine & Payoff Matrix

This section represents the computational implementation of Koulouris & Campajola’s memory hypothesis. The two agents recursively traverse 2 years of rolling historical windows (window_data).

At each node, they sample independent actions based on their weights, facing a non-cooperative game matrix:

• Mutual Cooperation (Action 0, 0): High joint payout (+10) mimicking a stable, supra-competitive margin.
• Mutual Aggressive Competition (Action Match): Low joint rent (+1), representing the competitive Nash baseline.
• Cheating / Under-cutting: Asymmetric penalization (+5 vs -5).

# 6. JOINT TRAINING ENGINE (Symmetric Memory Interaction)
message("Joint Training: Agent A & Agent B are learning Silver Market dynamics...")

for(i in 1:(nrow(silver_full) - T_horizon)) {
  window_data <- silver_full[i:(i + T_horizon - 1), ]
  vol <- sd(window_data$returns, na.rm = TRUE)
  if(is.na(vol)) vol <- 0.01
  
  state_vec <- matrix(c(1.0, vol, 0.5), nrow = 1)
  
  act_A <- which.max(predict(agent_A, state_vec, verbose = 0)) - 1
  act_B <- which.max(predict(agent_B, state_vec, verbose = 0)) - 1
  
  rewards <- if(act_A == 0 && act_B == 0) {
    list(A = 10, B = 10) 
  } else if(act_A == act_B) {
    list(A = 1, B = 1)   
  } else {
    if(act_A > act_B) list(A = 5, B = -5) else list(A = -5, B = 5) 
  }
  
  target_A <- predict(agent_A, state_vec, verbose = 0)
  target_B <- predict(agent_B, state_vec, verbose = 0)
  
  target_A[1, act_A + 1] <- rewards$A
  target_B[1, act_B + 1] <- rewards$B
  
  agent_A %>% fit(state_vec, target_A, epochs = 1, verbose = 0)
  agent_B %>% fit(state_vec, target_B, epochs = 1, verbose = 0)
}

7. Post-Convergence Audit Inference & Regime Selection

Once the networks stabilize, the engine takes the posture of an unbiased financial regulator. It extracts the neural policy configurations, evaluates the actual current execution window, and automatically determines the market regime using an automated classification layer.

# 7. FINAL AUDIT INFERENCE
analysis_data <- analysis_data %>%
  rowwise() %>%
  mutate(
    state_v = list(matrix(c(close/twap_path, current_sigma, (T_horizon - row_number())/T_horizon), nrow = 1)),
    q_A = list(predict(agent_A, state_v[], verbose = 0)),
    q_B = list(predict(agent_B, state_v[], verbose = 0)),
    joint_action = (which.max(q_A[]) + which.max(q_B[])) / 2
  ) %>% ungroup()

# 8. STATUS LOGIC (Professional Category Selection & Color Alignment)
last_row <- tail(analysis_data, 1)
market_status <- case_when(
  last_row$close >= last_row$twap_path ~ 
    list(
      label = "**COOPERATIVE:** Pareto-Efficient Alignment", 
      bg    = "#E8F8F5",  
      color = "#27AE60"   
    ),
  
  last_row$close < last_row$twap_path & last_row$close >= last_row$nash_path ~ 
    list(
      label = "**NORMAL:** Competitive Nash Equilibrium", 
      bg    = "#FEF5E7",  
      color = "#E67E22"   
    ),
  
  TRUE ~ 
    list(
      label = "**LIQUIDITY SHOCK:** Strategic Deviation Detected", 
      bg    = "#FDEDEC",  
      color = "#C0392B"   
    )
)

8. High-Fidelity Infographic Layer

To generate a publication-quality static vector infographic, we map our theme directly via ggplot2 and ggtext. By embedding the color palette directly into the HTML subtitle strings and forcing label formatting via scales::percent, we create a clean, high-contrast dashboard visualization.

# 9. GGPLOT PRODUCTION VISUALIZATION (Static Mode with ggtext Integration)
ggplot(analysis_data, aes(x = date)) +
  geom_ribbon(aes(ymin = lower_safety_limit, ymax = twap_path), fill = "darkgray", alpha = 0.3) +
  
  geom_line(aes(y = twap_path, color = "TWAP (Cooperative)"), size = 1) +
  geom_line(aes(y = nash_path, color = "Nash (Competitive)"), size = 1) +
  geom_line(aes(y = close, color = "Actual Price"), size = 1.3) +
  scale_y_continuous(labels = scales::label_currency()) +
  
  geom_richtext(
    aes(x = median(date), y = max(close, twap_path) * 1.02, label = market_status$label),
    fill = market_status$bg, color = market_status$color, size = 4,
    family = "Roboto Slab" 
  ) +
  
  scale_color_manual(
    name = NULL,
    values = c("Actual Price" = "steelblue", "TWAP (Cooperative)" = "#27AE60", "Nash (Competitive)" = "#E67E22")
  ) +
  
  labs(
    title = "Silver Market Strategic Audit Engine",
    subtitle = paste0(
      "<span style='color:#27AE60;'>─── **Cooperative Zone**</span> | ",
      "<span style='color:#E67E22;'>─── **Competitive Zone**</span> | ",
      "<span style='color:steelblue;'>─── **Actual Execution**</span><br><br>",
      "<span style='color:darkgrey;'>**Strategic Corridor** (Supra-Competitive Margin Zone)</span>"
    ),
    x = NULL, y = NULL,
    caption = glue("Dynamic Sigma: {scales::percent(current_sigma, accuracy = 0.01)} | Shortfall: {round(actual_cost, 2)}%")
  ) +
  
  theme_minimal(base_family = "Roboto Slab") +
  theme(plot.title = element_text(face = "bold", size = 16),
        plot.subtitle = element_markdown(face = "bold"), 
        axis.text = element_text(face = "bold"),
        legend.position = "none")

9. Empirics & Compliance Conclusion

When we run the complete inference loop on our terminal Silver execution window, the strategic narrative clarifies perfectly: Actual Execution (the blue trajectory) tracks downward, bypassing the cooperative upper envelope and adhering directly to the competitive boundaries.

The audit badge cleanly returns a status of NORMAL: Competitive Nash Equilibrium, with the terminal metrics computing the exact execution shortfall at 1.59% as indicated in the chart above. While the agents are technically complex neural networks capable of learning memory patterns, the actual price action during this specific ten-day horizon reflects a highly competitive regime, keeping the execution within standard Nash boundaries rather than shifting into a supra-competitive zone.

For quantitative auditors and systemic risk monitors, this approach signals a paradigm shift. Static threshold tests are blind to multi-agent learning trends. By deploying neural simulation baselines, structural compliance teams can automatically audit execution algorithms, isolating algorithmic alignment from pure market variance.

Green tests don’t automatically mean correct code.

Tests can be written in ways that pass reliably, survive every devtools::test() run, and still tell you almost nothing about whether your code works. The software engineering community has a name for these patterns: test smells. They’re structural problems in test code — not bugs, but patterns that signal the tests are doing less or something different than they appear.

Have you ever reviewed a pull request and only skimmed the test code, because why bother if you see in CI they’re green?

The problem with test code is that:

• it can be written in ways that it looks fine but is actually misleading,
• you may not review test code with the same seriousness as production code,
• you may not know what to look for in test code.

Here’s a list of problems in test code that show up regularly, and what to do about them.

1. Mystery Guest

🫣 The smell: the test relies on state that isn’t visible inside the test itself.

test_that("report uses the correct currency symbol", {
  result <- format_currency(1000)
  expect_equal(result, "£1,000")
})

This test passes on one machine and fails on another — and the reason isn’t in the test. Somewhere upstream, a .Rprofile file sets options(currency = "GBP"), or maybe it’s set with an environment variable. In a fresh session without that state, format_currency falls back to "$1,000".

The test has an invisible guest. You can’t understand it by reading it.

The fix is to make the dependency explicit — set it and clean it up within the test itself using withr:

test_that("report uses the correct currency symbol", {
  withr::with_options(list(currency = "GBP"), {
    result <- format_currency(1000)
    expect_equal(result, "£1,000")
  })
})

It might be also a default value in the function itself that is critical to the behavior being tested. In that case, the test should set it explicitly:

test_that("report uses the correct currency symbol", {
  # Arrange
  currency <- "GBP"
  amount <- 1000

  # Act
  result <- format_currency(amount, currency = currency)

  # Assert
  expect_equal(result, "£1,000")
})

Now the test tells its own story. Anyone reading it knows exactly what conditions it relies on.

2. Eager Test

🫣 The smell: one test that asserts on multiple unrelated behaviors.

test_that("data pipeline works", {
  result <- run_pipeline(raw_data)
  expect_equal(nrow(result), 100)
  expect_equal(ncol(result), 5)
  expect_true(all(!is.na(result$value)))
  expect_equal(result$category, expected_categories)
  expect_true(is.numeric(result$value))
})

This test is checking the row count, the column count, the NA handling, the category values, and the column type — all under the title “data pipeline works.” When it fails, the output tells you something broke. Not what.

Split it into cases where each title describes the specific behavior:

test_that("pipeline returns 100 rows", {
  # Arrange
  raw_data <- ...

  # Act
  result <- run_pipeline(raw_data)

  # Assert
  expect_equal(nrow(result), 100)
})

test_that("pipeline removes rows with missing values", {
  # Arrange
  raw_data <- ...

  # Act
  result <- run_pipeline(raw_data)

  # Assert
  expect_true(all(!is.na(result$value)))
})

test_that("pipeline produces numeric values", {
  # Arrange
  raw_data <- ...

  # Act
  result <- run_pipeline(raw_data)

  # Assert
  expect_true(is.numeric(result$value))
})

This test suite has more lines, yes. Is it bad? No.

We’re not optimising memory, we’re optimising for clarity. It’s also much easier to maintain — when a requirement changes, you know exactly which test to update and why.

Adding # Arrange, # Act, and # Assert comments is a simple habit that makes this smell hard to miss. When the # Assert block has five expect_* calls testing unrelated things, the label might prompt you to notice — and question whether they all belong together. When the # Act creates multiple results, it becomes clear that the test is doing too much.

The structure doesn’t fix the test, but it makes the problem easier to spot at a glance. The test is no longer a blob of code doing everything at once. It has a structure that invites you to think about what belongs in each section.

3. Over-specification

🫣 The smell: the test checks how the code does something, not what it produces.

test_that("sends notification on upload", {
  # Arrange
  mock_notify <- mockery::mock()
  mockery::stub(process_upload, "notify_user", mock_notify)

  # Act
  process_upload(file = "report.csv", user_id = "u42")

  # Assert
  mockery::expect_called(mock_notify, n = 1)
  mockery::expect_call(
    mock_notify,
    1,
    notify_user("u42", type = "upload_complete")
  )
})

This test will break the moment you rename notify_user, change its signature, or switch to a different notification mechanism internally — even if the behavior (the user gets notified) stays exactly the same.

Over-specified tests are tightly coupled to implementation. They punish refactoring.

The alternative is to test from the outside:

test_that("user receives notification after upload", {
  # Arrange
  notifications <- list()

  # Act
  process_upload(
    file = "report.csv",
    user_id = "u42",
    notify_user = function(user_id, ...) {
      notifications[[length(notifications) + 1]] <<- user_id
    }
  )

  # Assert
  expect_length(notifications, 1)
  expect_equal(notifications[[1]], "u42")
})

The test no longer cares how notification happens internally — just that "u42" received one. You can refactor the internals freely.

Notice what changed structurally: instead of patching notify_user after the fact with a stub, process_upload now accepts it as an argument. This is dependency injection, and it’s worth understanding as a pattern in its own right — not just as a way to avoid over-specification, but as a replacement for stubs and mocks altogether.

When a dependency is injected, the test passes in exactly the implementation it wants. There’s no patching, no intercepting, no targeting of internal function names. The test controls the dependency directly:

process_upload <- function(file, user_id, notify_user = real_notify_user) {
  # ... upload logic
  notify_user(user_id, type = "upload_complete")
}

In production, real_notify_user is used. In tests, any function with the same interface can be passed in. The key word is interface: the argument name and what it receives. That’s the contract between process_upload and whatever handles notifications.

This changes how test failures work — in a good way. With a stub, the test fails when the internal call site changes: rename the function, move it to a different module, inline it, and the stub can no longer target it. The failure is about test plumbing, not about behavior. With dependency injection, the test only fails when the interface changes — when notify_user is expected to receive different arguments, or when it’s removed entirely. That’s a failure you want to know about. It means the contract between components changed, and the fake dependency needs updating to match.

Notice also that the fake notify_user is just a plain function — no mocking library needed. A mock that accepts any call and records nothing will let over-specification creep back in through expect_called() and expect_call(). A simple fake that does one concrete thing — appends to a list, increments a counter, writes to a local variable — can only be asserted on from the outside. It doesn’t offer the machinery to inspect internals, so the test can’t become over-specified even if you wanted it to.

4. Sensitive Equality

🫣 The smell: the test asserts on an entire object when only part of it matters.

test_that("filters inactive users", {
  # Arrange
  users <- ...

  # Act
  result <- filter_active_users(users)

  # Assert
  expect_equal(result, data.frame(...))
})

The expectation was created when this test was first written. Since then, the users data frame has gained two new columns. Every test asserting against the full snapshot now fails — not because the filtering logic is wrong, but because the expected object is stale.

Asserting on a whole object isn’t always wrong.

A test that documents the full shape of the output — its columns, types, and structure — is genuinely useful. It acts as a contract: this is what filter_active_users returns. When that contract changes intentionally, the test tells you.

This matters especially in R. In a statically typed language the compiler enforces what a function returns — change the return type and the code won’t build. R has no such guarantee. A function can return a data frame today and a list tomorrow, drop a column, silently change a type, and nothing outside the tests will catch it.

A full-object contract test fills the loose-typing gap, giving you the kind of output verification the language doesn’t provide for free.

The smell appears when every test does this, or when the full-object assertion is used to check a behavior where most of the fields are irrelevant. If you’re testing that the filter removed inactive users, the output’s column names and data types are noise — and noise breaks tests for the wrong reasons.

A cleaner split is one test for the shape, separate tests for the behaviors:

test_that("returns a data frame with user_id and active columns", {
  # Arrange
  users <- ...

  # Act
  result <- filter_active_users(users)

  # Assert
  expect_equal(result, expected_users)
})

test_that("filters inactive users", {
  # Arrange
  users <- ...

  # Act
  result <- filter_active_users(users)

  # Assert
  expect_equal(result$user_id, c("u01", "u03", "u07"))
  expect_true(all(result$active))
})

The first test owns the full-object assertion and will break when the output shape changes — which is the right test to break. The second test owns the filtering logic and only breaks when the filtering logic is wrong – not when an irrelevant column is added or renamed.

5. Obscure Test

🫣 The smell: the setup is so long that the behavior being tested is buried inside it.

test_that("calculates discount correctly", {
  user <- list(
    id = "u42",
    name = "Alice",
    email = "alice@example.com",
    created_at = as.Date("2020-01-01"),
    plan = "pro",
    region = "EU",
    active = TRUE,
    referral_code = "REF001",
    payment_method = "card",
    currency = "EUR"
  )
  product <- list(
    id = "p99",
    name = "Annual Subscription",
    base_price = 200,
    category = "software",
    tax_rate = 0.2,
    available = TRUE
  )
  result <- calculate_discount(user, product)
  expect_equal(result$discount_pct, 0.15)
})

By the time you reach the assertion, you’ve waded through 20 lines of setup. Why does the user need a referral_code? Does currency affect the discount? Is tax_rate relevant? The intent is invisible.

The fix is to extract a builder that encodes sensible defaults and lets each test express only what’s relevant to that case:

# If only used in this test file, keep alongside the tests.
# If shared across multiple files, move to a `setup-*.R` file.
make_user <- function(...) {
  defaults <- list(
    id = "u_default", name = "Test User",
    plan = "basic", region = "EU", active = TRUE,
    currency = "EUR", payment_method = "card"
  )
  modifyList(defaults, list(...))
}

test_that("pro users receive a 15% discount", {
  # Arrange
  user <- make_user(plan = "pro")
  product <- list(
    id = "p99",
    base_price = 200,
    category = "software"
  )

  # Act
  result <- calculate_discount(user, product)

  # Assert
  expect_equal(result$discount_pct, 0.15)
})

Now the test title and the setup agree: this is a test about the "pro" plan. Everything else is noise that’s been moved out of sight.

Using # Arrange, # Act, and # Assert comments here also helps. The test is no longer a blob of code doing everything at once. It has a structure that invites you to think about what belongs in each section, and helps you understand when one secttion ends and another begins. The test is easier to read, and the intent is clearer.

6. Meaningless Title

🫣 The smell: the test name describes what the code does, not what behavior it verifies.

test_that("calculate_discount works", {
  # Arrange
  user <- make_user(plan = "pro")

  # Act
  result <- calculate_discount(user, product)

  # Assert
  expect_equal(result$discount_pct, 0.15)
})

This is by far the most common test smell I see.

Naming things is hard. I don’t blame anyone for struggling to find a good test title. But a meaningless title is as useless as no title at all. Let’s make some effort here! If we can pause to come up with a good name for a function — we can think of a name for a test too.

"calculate_discount works" is not a test title — it’s a function name with “works” appended. When this test fails, the output tells you that calculate_discount stopped working. It doesn’t tell you what it’s supposed to do or under what conditions it broke.

The title is the first thing you read when a test fails. It should answer two questions: what behavior is expected, and under what circumstances.

test_that("pro users receive a 15% discount", {
  # Arrange
  user <- make_user(plan = "pro")

  # Act
  result <- calculate_discount(user, product)

  # Assert
  expect_equal(result$discount_pct, 0.15)
})

Now the failure output reads: pro users receive a 15% discount — and you immediately know whether the logic is wrong or the test is outdated.

This smell is easy to spot: if your test title could be swapped with any other test on the same function without anyone noticing, it’s not a title — it’s a placeholder. A good title is a specification. It should still make sense even if the implementation is completely rewritten.

Some patterns that signal the smell:

• Ends with “works”, “is correct”, “should work”, or “test”
• Is just the function name
• Doesn’t mention the input condition that makes this case distinct

A rule of thumb: write the title before you write the test body. If you can’t name what you’re testing, you don’t know what you’re testing yet.

7. Flaky Test

🫣 The smell: the test sometimes passes and sometimes fails without any code change.

test_that("model selects the most important features", {
  # Arrange
  training_data <- ...
  top_n <- 3

  # Act
  result <- select_features(training_data, top_n = top_n)

  # Assert
  expect_equal(result$feature, c("age", "income", "region"))
})

This passes on the developer’s machine today. Tomorrow it returns c("income", "age", "region") — same features, different order — because two of them had near-identical importance scores and a random tie-breaking step inside the model produced a different ranking. The test fails. Nothing is broken.

Flaky tests are one of the most corrosive smells in a test suite. They don’t just fail — they train developers to distrust failures. Once the team learns to re-run tests until they go green, a genuinely broken build can hide for days.

The R-specific sources of flakiness to watch for: unset random seeds, tests with hard-coded Sys.sleep() waits that assume the app loads within a fixed time, tests that rely on network availability, and assertions against system time or floating-point values without a tolerance.

For the feature selection example, the fix is to set the random seed before the call, and to assert on what actually matters — that the right features were selected — rather than the incidental order they came back in:

test_that("model selects age, income, and region as top features", {
  # Arrange
  training_data <- ...
  top_n <- 3
  withr::local_seed(42)

  # Act
  result <- select_features(training_data, top_n = top_n)

  # Assert
  expect_equal(result$feature, c("age", "income", "region"))
})

Assuming it’s randomness that makes the test flaky, setting the seed makes any remaining randomness deterministic. The test now passes or fails for the right reason.

8. Never-Failing Test

🫣 The smell: the test will pass no matter what the code does.

test_that("returns an error for invalid input", {
  expect_error(calculate_discount(user = NULL, product))
})

This passes today. It will also pass if calculate_discount is rewritten to throw a completely different error, emit a warning instead, return NA silently, or crash with an unrelated message. The assertion only checks that some error occurred — not the right one.

The same trap appears with positive assertions:

test_that("result is a list", {
  result <- run_pipeline(data)
  expect_true(is.list(result))
})

If run_pipeline returns any list — empty, malformed, missing all expected fields — this test will almost always pass. It is not testing behavior; it is testing that R’s type system is still working.

Always assert on the specific outcome. For errors, check the message:

test_that("returns an error when user is NULL", {
  # Arrange
  user <- NULL
  product <- ...

  # Act, Assert
  expect_error(
    calculate_discount(user = NULL, product),
    "user must not be NULL"
  )
})

For positive assertions, check the values that actually matter — not just the type that contains them.

9. Test Order Dependency

🫣 The smell: tests only pass when run in a specific sequence.

test_that("initialises the cache", {
  initialise_cache()
  expect_true(cache_is_ready())
})

test_that("stores a value in the cache", {
  store_in_cache("key", "value")
  expect_equal(retrieve_from_cache("key"), "value")
})

The second test silently depends on the first having run. If you run it in isolation it fails with an error about the cache not being initialized, which looks like a cache bug but is actually a test structure bug.

The symptom is tests that pass with devtools::test() and fail with testthat::test_file() on a single file, or tests that break when you move them to a different describe() block.

The fix is to make each test fully self-contained. Every test should set up whatever it needs and leave no state behind:

test_that("stores a value in the cache", {
  # Arrange
  initialise_cache()
  withr::defer(clear_cache())
  store_in_cache("key", "value")

  # Act
  value <- retrieve_from_cache("key")

  # Assert
  expect_equal(value, "value")
})

withr::defer() ensures the cleanup runs even if the test fails, so the next test always starts from a clean slate.

10. Conditional Test Logic

🫣 The smell: the test body contains if, else, or tryCatch.

test_that("formats output correctly", {
  result <- format_output(data)
  if (nrow(result) > 0) {
    expect_equal(result$label[[1]], "Total")
  } else {
    expect_equal(nrow(result), 0)
  }
})

When a test branches, it is multiple tests in one — and only one branch executes on any given run. The other path is permanently untested. If the data changes so that nrow(result) is always 0, the label check never runs again and a regression in label formatting will go undetected forever.

Conditional logic in tests is also a sign the test doesn’t know what it expects. A test that hedges with if hasn’t decided what the correct behavior is.

If you find yourself using a tryCatch, then it seems like you don’t know whether the code should throw an error or not. Tests must know if an error is expected.

Split the conditions into separate, unconditional tests with setups that guarantee the state each one needs:

test_that("includes a Total label when data is non-empty", {
  # Arrange
  non_empty_data <- ...

  # Act
  result <- format_output(non_empty_data)

  # Assert
  expect_equal(result$label[[1]], "Total")
})

test_that("returns empty output when data is empty", {
  # Arrange
  empty_data <- ...

  # Act
  result <- format_output(empty_data)

  # Assert
  expect_equal(nrow(result), 0)
})

Each test now executes its assertion unconditionally, every time.

11. Test Logic in Production

🫣 The smell: production code contains special branches that only exist to make tests pass.

process_payment <- function(
  amount,
  user_id,
  env = Sys.getenv("APP_ENV")
) {
  if (env == "test") {
    return(list(
      status = "ok",
      transaction_id = "test-txn-001"
    ))
  }
  # ... real payment logic
}

The test environment gets a hardcoded success response. Production gets the real code. The two paths diverge immediately, so the tests are no longer testing what ships.

This smell often starts as a quick fix — “I’ll just skip the real API call in tests” — and ends up as a permanent fixture that no one dares remove. The production branch accumulates changes – the test branch stays frozen. At some point the tests are validating behavior that the production code no longer exhibits.

The correct fix is dependency injection: pass the thing you want to replace as an argument, and substitute it in tests:

process_payment <- function(
  amount,
  user_id,
  payment_gateway = real_gateway
) {
  result <- payment_gateway(amount, user_id)
  list(status = result$status, transaction_id = result$id)
}

test_that("returns ok status on successful payment", {
  # Arrange
  fake_gateway <- function(amount, user_id) {
    list(status = "ok", id = "txn-001")
  }
  amount <- 100
  user_id <- "u42"

  # Act
  result <- process_payment(
    amount,
    user_id,
    payment_gateway = fake_gateway
  )

  # Assert
  expect_equal(result$status, "ok")
})

Production code has no awareness of tests. Tests control the dependency. Both paths exercise the same logic.

The Pattern

These smells fall into two categories. Some make your tests untrustworthy: they pass for reasons unrelated to whether the code works. Others make your tests unreadable: they pass and fail correctly but no one can tell what they mean.

A test is documentation. It tells the next developer — probably future you — what the code is supposed to do and under what conditions. A smelly test fails at that job before it ever fails at catching bugs.

The most useful question to ask while writing a test isn’t “will this pass?” but “if this fails at 2am, will I understand what broke and why?”

If the answer is no, the test has a smell worth fixing.

Here’s a quick reference to keep nearby:

None of these require a major rewrite to fix. Each one has a mechanical correction: extract the state, split the block, inject the dependency, narrow the assertion, introduce a builder, rename the case. The value is in learning to see them first.

Once you can spot these smells, you can start to eliminate them, one by one — and build a better test suite that tells you the truth about your code.

Apply this

Reading about smells is easy. Spotting them in your own suite is the hard part. The fastest way to make them stick is to point an AI agent at your real test code, have it find the smells, then fix the worst few yourself — that’s where the learning happens.

Open a test file in your AI coding agent (Claude Code, Cursor, Copilot Chat) and paste this prompt:

You are a senior R engineer reviewing test code for "test smells" — structural patterns that let tests pass while hiding bugs or obscuring intent.

Scope: audit the test file(s) I've shared. Also read the production code they exercise — several smells can't be judged from the test alone. Use testthat and withr idioms in your fixes.

Detect these 11 smells. For each, here's the tell and the fix I prefer:

1. Mystery Guest — depends on options / env vars / .Rprofile set elsewhere → set and scope it inside the test with withr::with_options / local_envvar, or pass the value in as an argument.
2. Eager Test — one test_that asserts many unrelated behaviors → split to one behavior per test_that, each titled with the behavior it checks.
3. Over-specification — asserts on internal calls via mockery::stub / mock (expect_called, expect_call) → inject a fake as an argument and assert on the observable result, not the call.
4. Sensitive Equality — full-object expect_equal where only part matters → assert the specific fields. Keep ONE deliberate full-object "shape contract" test; that one is legitimate, not a smell.
5. Obscure Test — long setup buries the point → extract a make_*() builder with modifyList defaults, set only the field under test; move shared builders to setup-*.R.
6. Meaningless Title — title is "<fn> works" or just the function name → rewrite to state the expected behavior and the condition; it should still make sense if the implementation were rewritten.
7. Flaky Test — depends on seed / time / network / order, or asserts exact floats or incidental ordering → set withr::local_seed(), add a tolerance, assert on what matters not the order it came back in.
8. Never-Failing Test — expect_error() with no expected message, or asserts only a type (is.list, is.numeric) → assert the specific error message and the actual values.
9. Test Order Dependency — relies on state another test left behind → make it self-contained; set up what it needs and clean up with withr::defer().
10. Conditional Test Logic — if / else / tryCatch in the test body → split into unconditional tests whose setup guarantees each branch's state.
11. Test Logic in Production — `if (env == "test")`-style branches in the source → remove them and use dependency injection so tests and production run one path.

Rules:
- Only flag clear instances. Don't invent smells to look thorough, and don't flag the legitimate patterns above (injected fakes, a single shape-contract test, AAA comments).
- Quote the offending lines and cite file:line for every finding.
- Never change what the production code does — only its testability. If a fix needs a signature change (dependency injection), say so and describe the new interface.

Output:
1. A triage table: smell | file:line | severity (high/med/low) | one-line why.
2. Then fix the highest-severity findings as before/after code blocks — the smallest change that removes the smell. Stop after three and ask before continuing if more remain.
3. Tell me how to re-run just these tests to confirm they still pass.

Before you commit a test, run it past this checklist:

• The title states the behavior and the condition, and would still make sense if the implementation were rewritten.
• Everything the test depends on is visible inside the test (or its withr / setup-*.R scaffolding).
• It asserts on a specific outcome — not just a type, and not just “some error happened.”
• It passes or fails for exactly one reason, and that reason lives in the test, not in the run order.

Want a checklist for the whole suite, not just one file? The R testing roadmap turns these habits into a step-by-step path you can work through.

TL;DR, don’t let your friends use LLMs for finding useless code in a project! Using Jarl instead is cheaper, more reliable, and won’t kill any kitten.

We (Hannah and Maëlle) share an appreciation for the unglamorous maintenance work we call upkeep. So when Claude highlighted some dead code in the tune package to Hannah, it was worth a mention during a chat. This led to Maëlle recommending Etienne Bacher’s Jarl, a fast linter for R written in Rust that can, among other things, detect unused functions.

Motivated by this conversation among other things, Maëlle wrote a whole blog post about Jarl and other tooling for R based on tree-sitter, including the above meme.

So naturally, when you get called out in a custom meme, you listen. So here’s the long(er) form of “Use Jarl!”, from the following sweep through the parsnip package by Hannah, with commentary from Maëlle.

Many thanks to Etienne Bacher for reviewing this post!

How Hannah used Jarl

I have the Jarl extension installed for my Positron but here I used the command line interface (CLI) of Jarl. I started off with running this command in the terminal:

jarl check .

but that gave me too much output.

It did include a helpful hint though, something like:

── Summary ──────────────────────────────────────
Found 296 errors.
3 fixable with the `--fix` option.
More than 15 errors reported, use `--statistics` to get the count by rule.

So one more time, with the --statistics option, to get an overview of what Jarl all flagged. The output (here with made up numbers) looked something like

jarl check . --statistics
  132 [ ] implicit_assignment
  109 [ ] internal_function
   34 [ ] vector_logic
   13 [ ] unused_function
    7 [*] numeric_leading_zero
    1 [*] any_is_na

Rules with `[*]` have an automatic fix.

Now that gives me a better overall picture. Maëlle: I really enjoyed learning about this --statistics flag, and followed the same strategy when I ran Jarl on the igraph R package!

I could have run jarl check . --fix to fix everything with an automatic fix but I wanted to learn a bit more about the rules to see if I wanted to apply them. So I worked my way through one rule at a time, with one git commit at a time. You can select rules like so:

jarl check . --select any_is_na

This separation of changes by rule also made it easier to review changes: one topic per commit. Maëlle: in igraph, I created one PR per rule or group of rules.

I first went through the ones with automatic fixes, then the ones without. The result of the clean-up: https://github.com/tidymodels/parsnip/pull/1356 Maëlle: in igraph, I went through rules starting with the ones with the least violations, ending with the ones that had the most hits.

parsnip has been around the block a few times and had accumulated a few functions flagged by Jarl as unused. I spotted release_bullets() in the list which I know we want to keep.¹ We use it for parsnip in the release process but we don’t use it in parsnip. I wasn’t expecting Jarl to catch onto this difference so I decided to review the list of flagged functions separately and split that off into its own issue (and PR). I ended up keeping one other function that was part of a set of functions provided in a standalone file but I removed the other ones.

No kittens were hurt!

Why use Jarl for refactoring

Obviously, because we are telling you you should! More seriously: Jarl is fast and free!

And here are a few more reasons beyond our success stories :

Jarl will find small smelly things

For instance, the outer_negation rule helps improve readability: !all(x) is easier to understand than any(!x). In igraph for example, thanks to that rule we changed

any(!names(options) %in% names(defaults))

Any names of options not in names of defaults (and actually a tad worse since we do not use %notin% yet, new in R 4.6.0)

into

!all(names(options) %in% names(defaults))

Not all names of options in names of defaults.

See also the rule in action in parsnip.

Generally, Jarl will make the codebase less smelly.

Jarl will find unused functions

It’s very easy for a package to accumulate functions that are now useless, sort of historical artefacts. And as they’re unused you might not even stumble upon them by chance. Jarl will flag them for you!

If Jarl notices such a function, check it is actually unused. If the function is actually useful, you might want to add an exception comment or some rule-specific configuration.

Jarl will find duplicated function definitions

Also a problem that happens in a package that wasn’t born yesterday! In this case, you need to pick one definition. Example in parsnip. Deleting unused lines is very satisfying.

Jarl will detect unreachable code

For instance, code that comes after a return() or a stop() in a function. Sometimes, you can simply delete that code (again, a feel good move!). Other times, the check is an opportunity for targeted refactoring!

Jarl can check your testthat code

As a bonus, because these rules are turned off by default, you can run

jarl check . --select TESTTHAT

This will apply all rules from the testthat group that help you use more specific expectations. For instance, expect_equal(length(x), 2) should be expect_length(x, 2). All these rules come with automatic fixes!

Example in parsnip, example in igraph.

Should your package depend on dplyr, there’s a group of (currently two) rules for this as well!

Conclusion

Automatic tools are extremely useful for guiding upkeep work. Thinking “let me clean up this repo” can be daunting , tools like Jarl provide a roadmap! They let you fix obvious problems, and at the same time, they can make you read some dusty corners of your codebase. Hannah’s work on parsnip inspired Maëlle to run Jarl on igraph (influencer influenced back!), and she too can confirm it was both useful and satisfying. Thanks Etienne for that cool tool!

Note that you could run Jarl checks manually as we did, or you could tell a LLM to run Jarl. To quote Hannah’s colleague Emil Hvitfeldt, “When coding or using coding agents, one way to improve your workflow is by using CLI tools with a very clear focus.”. 🫡

What’s even better with Jarl is that it is actively developed, which you can confirm by skimming through the rich changelog!

The goodpractice package has been recommended by rOpenSci since it was first started just over 10 years ago by Gábor Csárdi. We used to ask our editors to manually run goodpractice on all packages submitted to software peer-review, and then to ask authors to fix any notable issues flagged by the package. It is now integrated within our own pkgcheck system, and used to automatically identify any goodpractice issues with all new submissions. The package changed maintainers several times before the previous maintainers gave the green light for us to take over maintenance of the package two years ago (28th May 2024).

We’re really pleased to share that we’ve recently rolled out a host of updates and extensions to the package. These make it both easier to use, and more powerful. This was a collaborative effort between new package author, Athanasia Mo Mowinckel, current maintainer, Mark Padgham, and the generative AI tool Claude. We describe the process at the end of this tech note.

Easier control of checks

The goodpractice package is a convenience wrapper around several other packages including rcmdcheck, lintr, cyclocomp, and desc, along with a few hand-coded checks within the package itself. The previous version of goodpractice had a total of 230 checks, generally prefixed with the name of the package which defined each check. Checks were identifiable with the single function, all_checks(). The only way to control which checks were run was to pass a checks parameter to the main function as a character vector of the names of checks you wanted to run. So, for example, if you wanted to skip the rcmdcheck checks, you had to do something like this:

mychecks <- grepv("^rcmdcheck", all_checks(), invert = TRUE)
gp(checks = mychecks)

That wasn’t easy! And honestly, most people just didn’t bother.

This update finally makes it easier to control checks, by defining them in “groups”, which generally correspond to the separate packages that run the actual checks. The groups can be seen with the new all_check_groups() function, and controlled with the checks_by_group() function. Checks can be controlled by simply naming the check groups you want to run, like this:

gp(checks = checks_by_group("description", "lintr"))

New check groups

The documentation for the all_check_groups() function shows all the new check groups added in this update, including:

Many of these checks are powered by treesitter, the extremely efficient syntax parsing library used by GitHub — which is genuinely neat to have working in R. The number of checks has been increased from 230 to 338, with all listed in the output of the all_checks() function.

Improved check reporting

The goodpractice package has always been designed for console output. This updated version now provides immediate and consistent detail on check progress while running. The cli package is also used to give consistently formatted output for all check groups.

The large increase in checks also means that directly printing goodpractice results may fill several screens of output. This update also adds a groups parameter to the print() method, to enable printing check results only for specified groups. This means you can run the main gp() function once, then step through each check group by printing results only for that group, fixing those, and then moving on to the next. We’ve found this turns what could feel like an overwhelming wall of output into a manageable to-do list. For example:

x <- gp()
print(x, "description")
#> Aww! Shining package! Keep up the priceless work!

print(x, "namespace")
#> ── It is good practice to ─────────────────────────────────────────────────────────────────────────────────────────────────────────────
#>
#> ✖ remove or use internal functions that are defined but never called. Dead code increases maintenance burden.
#>
#> R/utils.R:3
#> R/utils.R:73
#>
#> ✖ define exported (user-facing) functions before internal helper functions within each R source file.
#>
#> R/api.R:85
#>
#> ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Notes on the development process

From Mo (Athanasia)

I’ll be honest with you: a lot of this contribution was done with an AI assistant (Claude Code using Opus 4.5) at my elbow. I want to say something about that, because I think the honest version is more useful than the marketing version. The way Mark approached reviewing this work was very thoughtful and helpful.

Based on his feedback and how it was structured, I’ve learned a lot about what a good code review looks like, and it’s been absolutely wonderful. As a person who has mainly worked on their own code-bases alone, I have adopted a style of working that often creates diff monsters (hundreds of files changed), which can make reviewing very hard. Indeed, my first PR to this project was such a monster, and Mark pushed back asking for several smaller PRs he could actually manage to review. This is also where I learned about git worktrees (Thank you Maëlle!) and asking Claude to split the work up into several worktrees for easier review, was actually the first substantial change to how I now approach working with an AI assistant. Through this collaboration I learned a lot about splitting my work into better sizes and chunks.

The shape of it was something like this:

I’d locate an issue I thought was meaningful to tackle, and provide Claude Code with the issue context and possibly an idea of how to solve it. Once it had a solution, I’d review it, push back where it felt off, and we’d iterate. Mark would then review the PR and point out the things I’d missed. Three-way collaboration, more or less, with the AI doing the typing and me doing the deciding.

What worked well was scaffolding. Drafting a new check, wiring up the prep step that goes with it, generating the test cases for the edge paths I’d otherwise forget — that kind of work compresses really nicely with an AI. I would often ask it to generate some tests first, then create code that would pass the tests. This way we had a clear idea of what we wanted the new code to do, and then solve it.

What worked less well was anything involving judgement I hadn’t articulated yet. The first few times Claude opened a PR for me, it left the test-plan checkboxes unticked — even when the tests already passed. That’s not wrong exactly, but it’s misleading to a reviewer. I had to say “no, check the box if the test exists and passes; an unchecked box reads as TODO.” Similarly, when Mark left inline code suggestions on a PR, Claude tried to helpfully re-implement them locally — which would have stripped Mark’s attribution off his own contribution. We had to agree: suggestions get accepted on GitHub, not retyped.

And then there’s the unglamorous part. It would often forget instructions, despite having them documented in agents.md and local memory for the project. These were mostly trivial, but at times fairly bad. Some of the first PRs we made, Mark pointed out that the solutions were using regexp rather than AST (Abstrat Syntax Trees). I have to admit, I didn’t really know about AST before starting this project, and I am so glad I now know of it. Claude, despite being told very clearly that we wanted AST solutions, would often forget and implement regexp based solutions. It would also often revert to using for-loops rather than vectorization, and creating nested for-loops into the 3rd or 4th level — which is just horrible to follow as a human.

If you take one thing from this aside, take this: AI-assisted contribution didn’t mean handing the package over to a machine. It meant I could move faster on the parts I already understood, explore more confidently on the parts I didn’t, and spend my actual attention on the decisions that mattered — what to check for, how to group it, what to call it, and what to leave out. The friction taught me where to watch closely. The speed-ups taught me where I could trust the loop.

It’s a different way of working, and I’m still figuring out the shape of it, but it has enabled me to contribute to this project despite some very severe health issues making it hard for me to work as normal.

And one thing Claude always does better than me: writing good commit messages. I have found myself writing my own code, but asking Claude to commit them (making it look like Claude wrote the code, but I don’t care) because the commit messages are just so much better than what I write.

From Mark

The entire idea for this major update was initiated and driven by Mo. It was also my first real foray into using generative AI tools directly in rOpenSci packages, and I learnt a lot. After I rejected the initial monster PR, things settled down into much smaller and manageable PRs that were always focussed on one specific thing.

From that point, Mo was effectively the coder (assisted by Claude), and I was the reviewer. I never knew how much of the actual code was typed versus machine-generated, but was generally pleasantly surprised that I felt no need to inquire. Nearly every PR was obviously going to improve the package, and many of them were genuinely creative ideas. (That’s where I disagree slightly with Mo’s general description of the process: I am convinced a lot of the work behind this update came directly from her brilliant insights and experience, rather than her passive description of looking over pre-existing PRs.)

The process started in Feb 2026, and involved 70 pull requests. So many of these really were independent and new ideas. That enabled me to view each one with fresh eyes, and to think about whether I might approach anything in different ways. Over so many PRs, there was a lot of back-and-forth, from pushing back on Claude insisting on nested loops, to figuring out whether code was best parsed and analysed as text, or as a syntax tree. As Mo indicated, a lot of those discussions ultimately converged towards us bringing the power of treesitter into the package, and even helping to replace previous text-based code checks with more efficient and accurate AST approaches.

Our collaboration was very pleasant and enriching, and really felt like a direct cooperation between Mo as the ideas factory, me as the reviewer, and Claude as nothing more than a mediator of our ideas. The key role of Claude throughout the process was in handling all of the fiddly, technical details. We both looked carefully at every single line of code, but the use of Claude enabled our conversations to remain at higher, conceptual levels than what would have happened if we had to actually do all the technical implementation ourselves. I think that is the aspect that I found most surprising: That the use of Claude made our collaboration feel less technical, and therefore somehow even more human. And that gave us the ability to work though 70 pull requests representing over 100 new checks, all ready for everybody to use.

Let us know what you think

Like all rOpenSci packages, goodpractice is a community effort that lives through community use and feedback. We’d love to hear what you think, via email, or through issues on the GitHub repository.

https://www.techtonique.net contains tools for Exploratory Data Analysis (EDA), editors for R and Python code, tools for data visualization, no-code web interfaces for various data science tasks, a language-agnostic API for machine learning tasks (classification, regression, survival analysis, reserving, forecasting etc.).

I recently noticed a spike in sign-ups (not sure how they noticed it was up again, without official announcement…), and I thought it was a good opportunity to share that techtonique.net is back. More like a passion project, or part of a portfolio (for now). The API is available and evolving. You can use it for free (with rate limiting), but please note that it is slightly slower than it used to be (indeed, I used to have Azure Credits thanks to Microsoft for Startups => a faster server). It still serves the responses, in a reasonable time though, and I will try to optimize it as much as possible.

It is worth mentioning that I removed the stochastic simulation API for now, as it required (more technical) to run R in Python, through a Docker container.

Here are some examples of how to use the API for machine learning tasks in R and Python (there are also no-code interfaces for these tasks in the website, and the API is language-agnostic, so you can use it with any programming language that can make HTTP requests):

The starting point is to signup/login (if you’re facing issues, contact support@techtonique.net) and get a token from: https://www.techtonique.net/token. Then, you can use the token in API requests for machine learning tasks. Each response is a prediction, as envisaged by the chosen model, and based on the input data.

1 – R example:

#!/usr/bin/env Rscript

# Install httr if needed: install.packages("httr")
# All you need is a token from www.techtonique.net/token, then run the script
library(httr)

BASE_URL   <- "https://www.techtonique.net"
GITHUB_RAW <- "https://raw.githubusercontent.com/Techtonique/datasets/main"

DATASETS <- list(
  univariate     = paste0(GITHUB_RAW, "/time_series/univariate/a10.csv"),
  multivariate   = paste0(GITHUB_RAW, "/time_series/multivariate/ice_cream_vs_heater.csv"),
  classification = paste0(GITHUB_RAW, "/tabular/classification/breast_cancer_dataset2.csv"),
  regression     = paste0(GITHUB_RAW, "/tabular/regression/boston_dataset2.csv"),
  raa            = paste0(GITHUB_RAW, "/tabular/triangle/raa.csv"),
  abc            = paste0(GITHUB_RAW, "/tabular/triangle/abc.csv"),
  km             = paste0(GITHUB_RAW, "/tabular/survival/kidney.csv"),
  ridge_survival = paste0(GITHUB_RAW, "/tabular/survival/gbsg2_2.csv")
)

# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

get_token <- function() {
  args <- commandArgs(trailingOnly = TRUE)
  idx  <- which(args == "--token")
  if (length(idx) > 0 && idx < length(args)) return(args[idx + 1])
  cat("Please enter your JWT token: ")
  readLines(con = "stdin", n = 1)
}

# Download a GitHub dataset to a temp file and return its path.
# Using a real file is the only reliable way to do multipart uploads in httr.
fetch_dataset <- function(key) {
  url      <- DATASETS[[key]]
  filename <- basename(url)
  cat(sprintf("  Fetching %s from GitHub... ", filename))
  tmp <- tempfile(fileext = ".csv")
  resp <- GET(url, write_disk(tmp, overwrite = TRUE))
  stop_for_status(resp)
  cat("OK\n")
  tmp   # return the temp file path
}

make_request <- function(endpoint, token, dataset_key = NULL, params = list()) {
  url     <- paste0(BASE_URL, endpoint)
  headers <- add_headers(Authorization = paste("Bearer", token))

  if (!is.null(dataset_key)) {
    tmp_path <- fetch_dataset(dataset_key)
    on.exit(unlink(tmp_path), add = TRUE)   # clean up temp file when done
    body <- list(file = upload_file(tmp_path, type = "text/csv"))
    resp <- POST(url, headers, query = params, body = body, encode = "multipart")
  } else {
    resp <- GET(url, headers, query = params)
  }

  if (http_error(resp)) {
    cat(sprintf("  ERROR %s: %s\n", status_code(resp),
                content(resp, as = "text", encoding = "UTF-8")))
    return(NULL)
  }
  content(resp, as = "parsed", type = "application/json")
}

print_response <- function(resp) {
  if (is.null(resp)) return(invisible(NULL))
  for (key in names(resp)) {
    val <- resp[[key]]
    if (is.list(val) || length(val) > 6) {
      cat(sprintf("  %s: [%s ...]\n", key, paste(head(unlist(val), 6), collapse = ", ")))
    } else {
      cat(sprintf("  %s: %s\n", key, paste(val, collapse = ", ")))
    }
  }
}

# ---------------------------------------------------------------------------
# Test sections
# ---------------------------------------------------------------------------

test_forecasting <- function(token) {
  cat("\n=== Testing Forecasting Endpoints ===\n")

  cat("\nTesting univariate forecasting...\n")
  resp <- make_request("/forecasting", token,
    dataset_key = "univariate",
    params = list(base_model = "RidgeCV", n_hidden_features = 5,
                  lags = 25, type_pi = "kde", replications = 4, h = 3))
  print_response(resp)

  cat("\nTesting multivariate forecasting...\n")
  resp <- make_request("/forecasting", token,
    dataset_key = "multivariate",
    params = list(base_model = "RidgeCV", n_hidden_features = 5, lags = 25, h = 3))
  print_response(resp)
}

test_ml <- function(token) {
  cat("\n=== Testing Machine Learning Endpoints ===\n")

  cat("\nTesting classification...\n")
  resp <- make_request("/mlclassification", token,
    dataset_key = "classification",
    params = list(base_model = "RandomForestClassifier",
                  n_hidden_features = 5, predict_proba = TRUE))
  print_response(resp)

  cat("\nTesting regression...\n")
  resp <- make_request("/mlregression", token,
    dataset_key = "regression",
    params = list(base_model = "RidgeCV", n_hidden_features = 5, return_pi = TRUE))
  print_response(resp)
}

test_reserving <- function(token) {
  cat("\n=== Testing Reserving Endpoints ===\n")

  cat("\nTesting RidgeCV...\n")
  resp <- make_request("/mlreserving", token,
    dataset_key = "raa",
    params = list(method = "RidgeCV"))
  print_response(resp)

  cat("\nTesting LassoCV...\n")
  resp <- make_request("/mlreserving", token,
    dataset_key = "abc",
    params = list(method = "LassoCV"))
  print_response(resp)
}

test_survival <- function(token) {
  cat("\n=== Testing Survival Analysis Endpoints ===\n")

  cat("\nTesting Kaplan-Meier...\n")
  resp <- make_request("/survivalcurve", token,
    dataset_key = "km",
    params = list(method = "km"))
  print_response(resp)

  cat("\nTesting Ridge survival...\n")
  resp <- make_request("/survivalcurve", token,
    dataset_key = "ridge_survival",
    params = list(method = "RidgeCV", patient_id = 0))
  print_response(resp)
}

# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

main <- function() {
  token <- get_token()
  test_forecasting(token)
  test_ml(token)
  test_reserving(token)
  test_survival(token)
}

main()

2 - Python example:

import argparse
import requests
import os
import io
from typing import Dict, Any

# All you need is a token from www.techtonique.net/token, then run the script

# Base URL for API endpoints
base_url = "https://www.techtonique.net"

# Base URL for raw GitHub content
GITHUB_RAW = "https://raw.githubusercontent.com/Techtonique/datasets/main"

# Dataset paths within the GitHub repo
DATASETS = {
    "univariate":       f"{GITHUB_RAW}/time_series/univariate/a10.csv",
    "multivariate":     f"{GITHUB_RAW}/time_series/multivariate/ice_cream_vs_heater.csv",
    "classification":   f"{GITHUB_RAW}/tabular/classification/breast_cancer_dataset2.csv",
    "regression":       f"{GITHUB_RAW}/tabular/regression/boston_dataset2.csv",
    "raa":              f"{GITHUB_RAW}/tabular/triangle/raa.csv",
    "abc":              f"{GITHUB_RAW}/tabular/triangle/abc.csv",
    "km":               f"{GITHUB_RAW}/tabular/survival/kidney.csv",
    "ridge_survival":   f"{GITHUB_RAW}/tabular/survival/gbsg2_2.csv",
}


def get_token() -> str:
    """Get token from command line argument or prompt"""
    parser = argparse.ArgumentParser(description='Test API endpoints')
    parser.add_argument('--token', help='JWT token for authentication')
    args = parser.parse_args()

    if args.token:
        return args.token

    return input("Please enter your JWT token: ")


def fetch_dataset(dataset_key: str) -> tuple[io.BytesIO, str]:
    """
    Download a dataset from GitHub and return it as an in-memory bytes buffer
    together with a filename, ready to be passed to requests as a file upload.
    """
    url = DATASETS[dataset_key]
    filename = url.split("/")[-1]
    print(f"  Fetching {filename} from GitHub...", end=" ")
    response = requests.get(url)
    response.raise_for_status()
    print("OK")
    return io.BytesIO(response.content), filename


def make_request(
    url: str,
    token: str,
    method: str = "POST",
    file_tuple: tuple = None,   # (BytesIO, filename)
    params: Dict = None,
) -> Dict[str, Any]:
    """Make an API request and return the response."""
    headers = {"Authorization": f"Bearer {token}"}

    try:
        if method == "POST":
            files = None
            if file_tuple:
                buf, filename = file_tuple
                files = {"file": (filename, buf, "text/csv")}
            response = requests.post(url, headers=headers, files=files, params=params)
        else:
            response = requests.get(url, headers=headers, params=params)

        response.raise_for_status()
        return response.json()

    except requests.exceptions.RequestException as e:
        print(f"Error making request to {url}:")
        print(f"Status code: {e.response.status_code if hasattr(e, 'response') else 'N/A'}")
        print(f"Response: {e.response.text if hasattr(e, 'response') else str(e)}")
        return None


def test_forecasting(token: str):
    print("\n=== Testing Forecasting Endpoints ===")

    print("\nTesting univariate forecasting...")
    response = make_request(
        f"{base_url}/forecasting", token,
        file_tuple=fetch_dataset("univariate"),
        params={"base_model": "RidgeCV", "n_hidden_features": 5,
                "lags": 25, "type_pi": "kde", "replications": 4, "h": 3},
    )
    if response:
        print("Response:", response)

    print("\nTesting multivariate forecasting...")
    response = make_request(
        f"{base_url}/forecasting", token,
        file_tuple=fetch_dataset("multivariate"),
        params={"base_model": "RidgeCV", "n_hidden_features": 5, "lags": 25, "h": 3},
    )
    if response:
        print("Response:", response)


def test_ml(token: str):
    print("\n=== Testing Machine Learning Endpoints ===")

    print("\nTesting classification...")
    response = make_request(
        f"{base_url}/mlclassification", token,
        file_tuple=fetch_dataset("classification"),
        params={"base_model": "RandomForestClassifier",
                "n_hidden_features": 5, "predict_proba": True},
    )
    if response:
        print("Response:", response)

    print("\nTesting regression...")
    response = make_request(
        f"{base_url}/mlregression", token,
        file_tuple=fetch_dataset("regression"),
        params={"base_model": "RidgeCV", "n_hidden_features": 5, "return_pi": True},
    )
    if response:
        print("Response:", response)


def test_reserving(token: str):
    print("\n=== Testing Reserving Endpoints ===")

    print("\nTesting RidgeCV...")
    response = make_request(
        f"{base_url}/mlreserving", token,
        file_tuple=fetch_dataset("raa"),
        params={"method": "RidgeCV"},
    )
    if response:
        print("Response:", response)

    print("\nTesting LassoCV...")
    response = make_request(
        f"{base_url}/mlreserving", token,
        file_tuple=fetch_dataset("abc"),
        params={"method": "LassoCV"},
    )
    if response:
        print("Response:", response)


def test_survival(token: str):
    print("\n=== Testing Survival Analysis Endpoints ===")

    print("\nTesting Kaplan-Meier...")
    response = make_request(
        f"{base_url}/survivalcurve", token,
        file_tuple=fetch_dataset("km"),
        params={"method": "km"},
    )
    if response:
        print("Response:", response)

    print("\nTesting Ridge survival...")
    response = make_request(
        f"{base_url}/survivalcurve", token,
        file_tuple=fetch_dataset("ridge_survival"),
        params={"method": "RidgeCV", "patient_id": 0},
    )
    if response:
        print("Response:", response)


def main():
    token = get_token()
    test_forecasting(token)
    test_ml(token)
    test_reserving(token)
    test_survival(token)


if __name__ == "__main__":
    main()

Three hundred seventy-six of the new packages submitted to CRAN in April were still there in mid-May. Here are my Top 40 picks in twenty-three categories: Actuarial Analysis, Archaeology, Biology, Causal Inference, Computational Methods, Ecology, Economics, Environmental Studies, Epidemiology, Finance, Functional Data Analysis, Health Technology Assessment, Machine Learning, Medical Statistics, Meta Analysis, Networks, Physics, Programming, Statistics, Time Series, Utilities, and Visualization.

By default, ggplot2 places the legend outside the plot area (usually on the right). However, sometimes you may want to move the legend inside the plot to save space or improve the layout. This post explores how to achieve this using theme() parameters.

library(ggplot2)

ggplot(mtcars, aes(wt, mpg, color = factor(cyl))) +
  geom_point(size = 3) +
  theme_bw() +
  theme(legend.position = c(0.8, 0.8))

ggplot(mtcars, aes(wt, mpg, color = factor(cyl))) +
  geom_point(size = 3) +
  theme_bw() +
  theme(
    legend.position = c(1, 1),
    legend.justification = c("right", "top"),
    legend.background = element_blank()
  )

ggplot(mtcars, aes(wt, mpg, color = factor(cyl))) +
  geom_point(size = 3) +
  facet_wrap(~am) +
  theme_bw() +
  theme(
    legend.position = c(0.5, 0.2),
    legend.background = element_rect(fill = "white", color = "grey80")
  )

# Use 3 categories to leave the 4th spot in a 2x2 grid empty
mtcars_subset <- mtcars[mtcars$carb %in% c(1, 2, 4), ]

ggplot(mtcars_subset, aes(wt, mpg, color = factor(cyl))) +
  geom_point(size = 3) +
  facet_wrap(~carb, nrow = 2, ncol = 2) +
  theme_bw() +
  theme(
    # Coordinates c(0.75, 0.25) target the center of the empty bottom-right area
    legend.position = c(0.75, 0.25),
    legend.background = element_blank(),
    legend.justification = "center"
  )

© 2025 tidyomics. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

You’ve set up all the preconditions.

You’ve triggered the action.

Now the specification needs to answer one question: did the system do the right thing?

That’s what Then steps are for. How precisely you answer it determines whether your specifications act as a reliable safety net or produce false confidence.

This article is the 3rd part of a series on writing BDD specifications for Shiny applications. We’ve built a data submission form, managed preconditions with Given steps, and modeled user interactions with When steps.

Read the previous articles to get full context, or continue here to focus on writing Then steps.

1. Behavior-Driven Development in R Shiny: A Step-By-Step Example
2. Behavior-Driven Development in R Shiny: Setting Up Test Preconditions with Given Steps
3. Behavior-Driven Development in R Shiny: Writing When Steps That Model User Behavior

Level-up your testing game! Grab your copy of the R testing roadmap.

The Purpose of Then

Then steps answer one question: What changed as a result of the user’s action?

Not how it changed internally. What the user can now observe — or what the system has now done — that wasn’t true before.

This distinction shapes every Then step you write. The assertion lives at the level of the outcome, not at the level of the mechanism that produced it.

# Outcome — what the user cares about
then_there_are_entries <- function(context, n) {
  expect_equal(context$storage$size(), n)
  context
}

# Mechanism — how the app stored it internally
then_the_database_table_has_n_rows <- function(context, n) {
  expect_equal(nrow(DBI::dbReadTable(context$conn, "entries")), n)
  context
}

Both assertions verify the same underlying fact — that the entry was saved — but the first one survives a switch from disk cache to database to API. The second breaks the moment you change the storage technology.

Assert what the user or system guarantees. Leave how it’s achieved to implementation tests.

What to Assert

Then steps verify outcomes across three kinds of observable evidence. These categories describe what you can check — they don’t prescribe how many steps to use. How to group assertions into steps is a separate question, covered in the next section.

User-visible state

Did the UI change the way a user would expect?

then_i_am_prompted_to_provide_required_fields <- function(context) {
  context$driver$expect_validation_feedback()
  context
}

Assertions use driver methods — the same translation layer that When steps use — so they’re protected from UI implementation details.

System state

Sometimes the most meaningful evidence isn’t what the user sees but what the system now holds.

then_the_entry_has_title <- function(context, expected_title) {
  entry <- context$storage$get_first()
  expect_equal(entry$title, expected_title)
  context
}

These assertions reach directly into the storage object that Given steps configured. They’re faster than UI assertions and more precise: they check the exact state of the system without waiting for the browser to render anything.

Side effects

Some behaviors produce no visible UI change and no stored data — they trigger side effects: emails, API calls, log entries. These still need verification.

then_email_notification_is_sent <- function(context) {
  context$email_service$expect_sent()
  context
}

then_i_am_informed_email_was_not_sent <- function(context) {
  context$driver$expect_visible("email_failure_message")
  context
}

The first assertion works because the email service test double (set up in Given) records whether it was called. The second checks that the app communicated the failure to the user — which is the behavior the specification is actually about.

Implementing Then Steps in the Driver

Then steps that check UI state belong in the driver, just like When steps. This keeps the assertion details hidden from the specification:

#' tests/testthat/setup-driver.R
MyAppDriver <- R6::R6Class(
  classname = "MyAppDriver",
  inherit = shinytest2::AppDriver,
  public = list(
    expect_visible = function(output_id) {
      # Find if HTML element is visible
      invisible(self)
    },
    expect_entry_count = function(n) {
      # Find number of HTML elements
      invisible(self)
    },
    expect_validation_feedback = function() {
      # Check if HTML element is visible
      invisible(self)
    }
  )
)

Each driver method is a named, reusable assertion. If the HTML structure changes, you fix it in one place.

Grouping Then Steps

The right question is not “how many assertions per step?” but “which assertions belong together?”

Group assertions that describe the same observable behavior. Split when outcomes can diverge independently.

Consider what happens when a user successfully submits an entry: the storage grows by one, and a confirmation message appears. These two outcomes are two sides of the same coin — if one is true, the other should be too. Splitting them into separate steps implies they can diverge, which invites the reader to wonder what it would mean for storage to succeed but no confirmation to appear, or vice versa. Grouping them into a single step names the behavior directly:

then_entry_is_submitted <- function(context) {
  expect_equal(context$storage$size(), 1)
  context$driver$expect_visible("confirmation_message")
  context
}

Email notification is different. Submission can succeed even when the email service fails — the two outcomes genuinely can diverge. That’s exactly when a separate step is the right call:

it("should submit entry and send notification", {
  given_no_content() |>
    given_an_authenticated_user(email = "user@example.com") |>
    given_email_service_is_available() |>
    when_i_submit_entry_with_all_required_fields() |>
    then_entry_is_submitted() |>
    then_notification_was_sent_to_the_authenticated_user()
})

Two steps, not three — because the spec now reflects the actual structure of the behavior. If then_entry_is_submitted fails, you know the core submission broke. If then_notification_was_sent_to_the_authenticated_user fails, you know the side effect broke. The split carries information because it maps to a real divergence point.

Then steps should also be read-only. They inspect state; they don’t change it. A Then step that modifies storage or triggers side effects is doing the wrong job. Keep the flow clean: Given sets up, When acts, Then observes.

Making Failure Messages Helpful

A failing test with a good message saves minutes. A failing test with a bad message wastes hours.

The most common mistake is letting low-level assertion failures surface directly. When expect_equal(nrow(df), 2) fails with "actual 0, expected 2", that tells you nothing about which scenario failed or what the data looked like.

The label argument sets the name testthat uses for the object (first argument) in the failure message. Keep it short and descriptive — testthat appends the actual and expected values itself:

then_there_are_entries <- function(context, n) {
  testthat::expect_equal(
    context$storage$size(), n,
    label = "number of entries in storage"
  )
  context
}

A failure now reads: "Expected number of entries in storage to equal n" — immediately clear without any manual formatting.

For UI assertions, label the element being checked:

expect_visible = function(output_id) {
  val <- self$get_value(output = output_id)
  testthat::expect_true(
    !is.null(val) && nchar(val) > 0,
    label = sprintf("output '%s'", output_id)
  )
  invisible(self)
}

A failure reads: "Expected output 'confirmation_message' to be TRUE" — which immediately tells you which element to look at.

What Doesn’t Belong in Then Steps

Don’t push implementation tests up to Then steps.

If you’re asserting the exact text of an error message, the CSS class applied to an invalid input, or the exact SQL that was executed — those are implementation details. They don’t belong here.

Acceptance-level Then steps answer: did the right thing happen from the user’s perspective?

The Complete Picture

Putting it all together, a well-formed specification has a clear structure at every level:

describe("data submission", {
  it("should submit entry and send notification", {
    given_no_content() |>
      given_an_authenticated_user(email = "user@example.com") |>
      given_email_service_is_available() |>
      when_i_submit_entry_with_all_required_fields() |>
      then_entry_is_submitted() |>
      then_notification_was_sent_to_the_authenticated_user()
  })

  it("should handle email service failure gracefully", {
    given_no_content() |>
      given_an_authenticated_user() |>
      given_email_service_is_unavailable() |>
      when_i_submit_entry_with_all_required_fields() |>
      then_entry_is_submitted() |>
      then_i_am_informed_email_was_not_sent()
  })

  it("should require all required fields", {
    given_no_content() |>
      given_an_authenticated_user() |>
      when_i_submit_entry_with_missing_required_fields() |>
      then_i_am_prompted_to_provide_required_fields()
  })
})

Each specification reads like a sentence. None of them mention shinytest2, input IDs, or database queries. The implementation lives behind the DSL, where it can change without touching the specifications.

Wrapping up

Then steps are where a specification earns its credibility.

A specification that only tests easy outcomes produces false confidence. One that asserts user-visible state, system state, and side effects actually catches real problems. Assert outcomes, not implementations. Group assertions by observable behavior and split them only where they can genuinely diverge. Write failure messages that save time. Push implementation details down to unit and module tests.

With Given, When, and Then steps in place, the specifications read like requirements and run on every build. That’s worth more than any testing framework on its own.

I’ve been using Unix and Linux systems for nearly 20 years. When I found myself with a Windows laptop for a while I felt lost and quickly installed Windows Subsystem for Linux. And yet the number of commands I know and use regularly is surprisingly short: cd, ls, rm, pwd, mkdir, source, node, npm and, of course, git. Even then I don’t think I really “know” git.

In his excellent book “Don’t Make Me Think Revisited”, Steve Krug describes how

Faced with any sort of technology, very few people take the time to read instructions. Instead, we forge ahead and muddle through, making up our own vaguely plausible stories about what we’re doing and why it works.

Krug’s book is about web and mobile usability but the point he makes here applies across technological domains. And it is safe to say I’ve been “muddling through” with Unix. If you asked me whether something was a Unix command or part of the shell I’d probably just shrug and say something like “I… I don’t know that!?”.

Until recently StackOverflow had always been my crutch to fall back on. From the frequency with which it appears in answers, I know that find is a very useful and versatile command. However, the odds of me remembering how to use it without looking it up are indistinguishable from 0.

But now we have alternatives in the form of various large language models or, to be more precise, multi-modal AI Assistants (like Claude and ChatGPT) and Agents (like Claude Code and Codex). I’d describe myself as AI-agnostic: I see benefits, I worry about issues. So it was with some (I’d say healthy) degree of scepticism I started using Claude to assist me in writing some shell functions to automate tasks I do or would like to do regularly.

I say “I” wrote. Reality is of course that I asked Claude to write and rewrite things for me and then I questioned why it did what it did. This probably makes it sound like I was assessing Claude for a job but it’s actually that I want to understand what has been written and, as we’ve already covered, I am not very good at this stuff.

Task One: fixing a common git mistake

One of my most common slip-ups with git is to commit changes to one branch (usually main) that were meant to go on a new feature branch. I know a way to recover from this: create the new branch I meant to create all along; do an interactive rebase to remove the commit from main; switch to the new feature branch I just created. There’s not that many steps to the process but it’s still a bit of a faff. Let’s ask Claude:

# Create and switch to the correct branch from current state
git switch -c correct-branch

# Go back to the wrong branch and undo the commit (keeps changes staged)
git switch wrong-branch
git reset --soft HEAD~1

# Discard the changes
git restore --staged .
git restore .

git switch -c correct-branch
git switch wrong-branch
git reset --hard HEAD~1

Here we hit a couple of issues.

1. Neither of Claude’s suggestions did the final step of what I wanted: moving me back to the new, correct, branch. This is entirely on me. I forgot to ask (and didn’t notice for a bit). I do that kind of thing quite a lot.
2. Even for simple tasks, Claude’s responses can be quite verbose.

The problem here is not a practical one regarding the use of Claude or an alternative agent, the real problem is it’s quite difficult to write a blog post of a sensible length that actually goes anywhere while accurately recounting the conversation you have with Claude. So what follows is a paraphrase rather than a transcript; shortened for the sake of your sanity and mine.

Tim

Can I use git switch – instead?

Claude

Yes, and it’s cleaner. git switch - saves you typing the wrong branch name out explicitly.

Tim

Seems like the sort of thing I should make into a Bash function

Claude

Agreed. Something like this (with an added safety check):

git-move-last() {
 if [ -z "$1" ]; then
 echo "Usage: git-move-last <new-branch-name>"
 return 1
 fi
 git switch -c "$1"
 git switch -
 git reset --hard HEAD~1
}

Tim

I’ve gone for git-behead. If I wanted to move to the new branch at the end I just add another git switch -?

Claude

Ha .

Exactly:

git-behead() {
 if [ -z "$1" ]; then
 echo "Usage: git-behead <new-branch-name>"
 return 1
 fi
 git switch -c "$1"
 git switch -
 git reset --hard HEAD~1
 git switch -
}

Tim

Ok, one more idea. Suppose I haven’t been looking what I’m doing and I want to remove more than just the last commit to a new branch. Can we recreate the above but with an optional commit count that defaults to 1?

Claude

git-behead() {
 if [ -z "$1" ]; then
 echo "Usage: git-behead <new-branch-name> [commit-count]"
 return 1
 fi
 local count="${2:-1}"
 git switch -c "$1"
 git switch -
 git reset --hard HEAD~"$count"
 git switch -
}

Tim

Could you comment anything you think isn’t self-explanatory?

Claude

git-behead() {
 if [ -z "$1" ]; then
 echo "Usage: git-behead <new-branch-name> [commit-count]"
 return 1
 fi
 local count="${2:-1}" # default to 1 commit if not specified
 git switch -c "$1" # create new branch at current HEAD
 git switch - # switch back to original branch
 git reset --hard HEAD~"$count" # remove commits from original branch
 git switch - # switch to new branch
}

Putting this in my .bashrc file gives me a simple command that I can invoke whenever I make the mistake I regularly make and an optional argument to use whenever I haven’t been paying attention to the branch I’m on for an extended period. It’s also short enough that I think I follow every line, despite my limited knowledge of Bash, and there are comments there for when future me has forgotten. This conversation took about fifteen minutes and included more questions than I’ve included here, plus time I spent running various ad-hoc tests with the different iterations of the function. That seems like a pretty efficient use of my time.

Task Two: listing image files

I do a lot of data visualisation, I look after the artwork for this website and I take photographs for our conferences. In short, I deal with a lot of images in various formats. It would be nice to be able to list only images and their key properties in a specific location with a simple function. Let’s ask Claude:

Tim

Can you write me a Bash script for displaying a list of images in a directory that shows file name, width, height and size (in kB) and nothing else?

Claude

Here’s a clean script using identify from ImageMagick (widely available) with a du fallback for file size:

#!/usr/bin/env bash

DIR="${1:-.}"

if ! command -v identify &>/dev/null; then
 echo "Error: imagemagick is required (install with: brew install imagemagick / apt install imagemagick)" >&2
 exit 1
fi

shopt -s nullglob nocaseglob
files=("$DIR"/*.{jpg,jpeg,png,gif,webp,bmp,tiff,tif,avif})
shopt -u nullglob nocaseglob

if [[ ${#files[@]} -eq 0 ]]; then
 echo "No images found in: $DIR" >&2
 exit 0
fi

printf "%-40s %8s %8s %10s\n" "FILENAME" "WIDTH" "HEIGHT" "SIZE (kB)"
printf "%-40s %8s %8s %10s\n" "--------" "-----" "------" "---------"

for f in "${files[@]}"; do
 name="$(basename "$f")"
 read -r width height < <(identify -format "%w %h" "$f" 2>/dev/null | head -1)
 size_kb=$(awk "BEGIN { printf \"%.1f\", $(wc -c < "$f") / 1024 }")
 printf "%-40s %8s %8s %10s\n" "$name" "$width" "$height" "$size_kb"
done

Supported formats: jpg, jpeg, png, gif, webp, bmp, tiff, avif.

Example output:

FILENAME WIDTH HEIGHT SIZE (kB)
-------- ----- ------ ---------
banner.png 1200 630 84.3
avatar.jpg 400 400 32.1
hero.webp 1920 1080 210.7

I asked for a script, I got a script. I wanted a function. This one is, again, on me.

Tim

Sorry, can we make this a function rather than a script?

Claude

The key changes: shebang removed, exit replaced with return, and variables are now local. To use it, source it in your shell and call it as normal:

imgls() {
 local dir="${1:-.}"

 if ! command -v identify &>/dev/null; then
 echo "Error: imagemagick is required (install with: brew install imagemagick / apt install imagemagick)" >&2
 return 1
 fi

 shopt -s nullglob nocaseglob
 local files=("$dir"/*.{jpg,jpeg,png,gif,webp,bmp,tiff,tif,avif})
 shopt -u nullglob nocaseglob

 if [[ ${#files[@]} -eq 0 ]]; then
 echo "No images found in: $dir" >&2
 return 0
 fi

 printf "%-40s %8s %8s %10s\n" "FILENAME" "WIDTH" "HEIGHT" "SIZE (kB)"
 printf "%-40s %8s %8s %10s\n" "--------" "-----" "------" "---------"

 local f name width height size_kb
 for f in "${files[@]}"; do
 name="$(basename "$f")"
 read -r width height < <(identify -format "%w %h" "$f" 2>/dev/null | head -1)
 size_kb=$(awk "BEGIN { printf \"%.1f\", $(wc -c < "$f") / 1024 }")
 printf "%-40s %8s %8s %10s\n" "$name" "$width" "$height" "$size_kb"
 done
}

The conversation goes on a bit from here and the function becomes more and more detailed. It would be tedious to show the whole conversation but to summarise, “we”:

1. Add support for SVG, including extracting dimensions from the viewBox attribute.
2. Fix a bug with GIF dimensions.
3. Tweak the header text.
4. Add support for recursive searching of a directory.
5. Fix the layout when doing a recursive search.
6. Round the values in the SIZE column
7. Add comments.

The final result:

imgls() {
 local recursive=0

 # Parse flags — currently only -r (recursive) is supported
 while getopts ":r" opt; do
 case $opt in
 r) recursive=1 ;;
 *) echo "Usage: imgls [-r] [directory]" >&2; return 1 ;;
 esac
 done
 # Shift past any parsed flags so $1 is the directory argument
 shift $((OPTIND - 1))
 OPTIND=1

 # Use first argument as directory, default to current directory if omitted
 local dir="${1:-.}"

 # identify is provided by imagemagick and used to read raster image dimensions
 if ! command -v identify &>/dev/null; then
 echo "Error: imagemagick is required (install with: brew install imagemagick / apt install imagemagick)" >&2
 return 1
 fi

 local exts=("jpg" "jpeg" "png" "gif" "webp" "bmp" "tiff" "tif" "avif" "svg")
 local files=()

 if [[ $recursive -eq 1 ]]; then
 # Build a find command matching all extensions, case-insensitively
 local find_args=("$dir" "-type" "f" "(")
 for i in "${!exts[@]}"; do
 [[ $i -gt 0 ]] && find_args+=("-o")
 find_args+=("-iname" "*.${exts[$i]}")
 done
 find_args+=(")")
 while IFS= read -r f; do
 files+=("$f")
 done < <(find "${find_args[@]}" | sort)
 else
 # nullglob: unmatched globs expand to nothing instead of being treated as literals
 # nocaseglob: match extensions case-insensitively (e.g. .JPG as well as .jpg)
 shopt -s nullglob nocaseglob
 files=("$dir"/*.{jpg,jpeg,png,gif,webp,bmp,tiff,tif,avif,svg})
 shopt -u nullglob nocaseglob
 fi

 if [[ ${#files[@]} -eq 0 ]]; then
 echo "No images found in: $dir" >&2
 return 0
 fi

 # In recursive mode show the relative path, otherwise just the filename
 printf "%-40s %11s %12s %10s\n" "FILENAME" "WIDTH (px)" "HEIGHT (px)" "SIZE (kB)"
 printf "%-40s %11s %12s %10s\n" "--------" "----------" "-----------" "---------"

 local f name width height size_kb viewbox current_subdir="" subdir
 for f in "${files[@]}"; do
 if [[ $recursive -eq 1 ]]; then
 # Print a header line when we enter a new subdirectory
 subdir="$(dirname "${f#$dir/}")"
 if [[ "$subdir" != "$current_subdir" ]]; then
 [[ -n "$current_subdir" ]] && echo ""
 echo "${subdir}/"
 current_subdir="$subdir"
 fi
 name="$(basename "$f")"
 else
 name="$(basename "$f")"
 fi
 # wc -c counts raw bytes; awk converts to kB with one decimal place
 size_kb=$(awk "BEGIN { printf \"%d\", $(wc -c < "$f") / 1024 }")

 if [[ "${f,,}" == *.svg ]]; then
 # SVGs are XML so we can parse dimensions directly without imagemagick.
 # Prefer viewBox (e.g. "0 0 100 200") as it's the most reliable source
 # of the intended coordinate space; fall back to width/height attributes.
 viewbox=$(grep -oi 'viewBox="[^"]*"' "$f" 2>/dev/null | head -1 | grep -oi '"[^"]*"' | tr -d '"')
 if [[ -n "$viewbox" ]]; then
 # viewBox is "min-x min-y width height" — fields 3 and 4
 width=$(echo "$viewbox" | awk '{printf "%g", $3}')
 height=$(echo "$viewbox" | awk '{printf "%g", $4}')
 else
 width=$(grep -oi 'width="[^"]*"' "$f" 2>/dev/null | head -1 | grep -oi '"[^"]*"' | tr -d '"px"')
 height=$(grep -oi 'height="[^"]*"' "$f" 2>/dev/null | head -1 | grep -oi '"[^"]*"' | tr -d '"px"')
 fi
 # ${var:--} expands to "-" if var is empty or unset
 width="${width:--}"
 height="${height:--}"
 else
 # [0] selects only the first frame, preventing animated GIFs from
 # producing one line of output per frame
 read -r width height < <(identify -format "%w %h" "${f}[0]" 2>/dev/null)
 fi

 printf "%-40s %11s %12s %10s\n" "$name" "$width" "$height" "$size_kb"
 done
}

The “development time” here was longer and I’d be fibbing if I said I understood every line of this function. But I do get the general principles and understand what it can and can’t do. And being able to list all images (of nearly all relevant types) in a directory with just their key properties in a nice tabular format is a very nice productivity win for me.

Final Thoughts

Discussions about AI easily tend towards the extremes: ‘it’s amazing’, ‘it’s a disaster’, ‘it’s useless’. This is particularly unhelpful when focusing on how successfully or otherwise it executes a single task. All we’ve really done in this post is create a couple of Bash scripts. For me the AI in question, Claude, has been useful. No more, no less. Despite this, the process did make me think a lot about how the way we interact with machines is changing. While it’s easy to over-personify AI, this line from the film that gave this blog post its title seems now, in 2026, to have been prescient:

The advent of computers, and the subsequent accumulation of incalculable data has given rise to a new system of memory and thought parallel to your own.

For updates and revisions to this article, see the original post

To leave a comment for the author, please follow the link and comment on their blog: The Jumping Rivers Blog.

---

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Ghost in the Shell Script]]> 401507 https://www.r-bloggers.com/2026/05/ropensci-news-digest-may-2026/ Thu, 28 May 2026 00:00:00 +0000 https://ropensci.org/blog/2026/05/28/news-may-2026/

Dear rOpenSci friends, it’s time for our monthly news roundup! You can read this post on our blog. Now let’s dive into the activity at and around rOpenSci!

rOpenSci HQ

15 Years of rOpenSci, and we’re just getting started 🎉
This y...

Continue reading: rOpenSci News Digest, May 2026]]> [social4i size="small" align="align-left"]

-->

[This article was first published on rOpenSci - open tools for open science, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Dear rOpenSci friends, it’s time for our monthly news roundup! You can read this post on our blog. Now let’s dive into the activity at and around rOpenSci!

rOpenSci HQ

15 Years of rOpenSci, and we’re just getting started

This year we celebrate 15 years of rOpenSci! First up we have a coworking session followed by a couple casual virtual community celebrations. Please join us! Read more in our blog post about these upcoming events.

Also exciting is that this blog post about our 15 year anniversary was featured in the R Weekly Highlights podcast hosted by Eric Nantz and Mike Thomas.

Champions Program update

The rOpenSci Champions Program is currently running two very active cohorts, with many activities happening across the program, and everything is on track!

The 2025–2026 cohort is wrapping up their projects and outreach activities, with several Champions presenting their work at events such as R/Pharma, useR!, Posit, and LatinR 2026. We also recently held a cohort meeting bringing together Champions and mentors from the 2025–2026 cohort and one meeting for mentors only, to share mentoring tips and tricks.

Meanwhile, the 2026–2027 cohort has started their training activities, including three weekly training workshops for Champions and a mentor training workshop. Mentors and Champions have also had their first meetings together to get to know each other and begin planning their collaboration.

Stay tuned to our events page for Champions outreach activities and to our blog to learn more about Champions’ projects.

May is Maintainer Month

Open source software doesn’t sustain itself. Behind every R package, there is at least one person who responds to issues, reviews pull requests, keeps up with dependency changes, and makes sure everything still works.

We joined the #MaintainerMonth celebration with a month-long series of rOpenSci’s maintainer spotlights: 36 maintainers from 15 countries across 4 continents, maintaining more than 50 packages that together serve thousands of researchers and data practitioners around the world.

Follow our social media (Mastodon and LinkedIn) and the official hashtag (#MaintainerMonth) to learn more about the maintainers who support open source software. Thank you to those who agreed to participate and to all the maintainers in our community.

Software-Review update

This month saw a slight policy shift in the way editors handle new software review submissions, assisted by an internal update to our bot system. Previously, our rotating Editor-in-Chief (EiC) was solely responsible for initial decisions on whether to proceed with submissions. Once a decision was reached, a handling editor took over, and generally guided each submission through to approval. However, recent increases in submissions have been too much for one EiC to handle. Submissions will now be passed on more quickly to handling editors who will then make initial scope and fit decisions.

As part of this change, the EiC can now ask the bot to “ping editors”. Editors will receive an email alerting them of this new submission and can then decide whether to take it on as editor. Each time an editor clicks on a link in the email, the EiC will receive an email notification, and once an editor is assigned, all links will be automatically deactivated.

For anybody interested in how to develop an email alert system like this, please contact us and we’ll have happy to share what we learnt along the way. Full code is in this R file.

Social media cards for R-universe

When posting about your package or universe on social media, you can get a pretty social media card featuring important information about (respectively):

• the package such as its maintainer, its keywords, its URL, its version, its number of stars, etc.

• the universe such as its number of packages and contributors, etc.

More information in the R-universe documentation.

Coworking

Read all about coworking!

• Tuesday June 2nd 2026, 14:00 Europe Central (12:00 UTC) “15 years with rOpenSci”, with Steffi LaZerte and cohost Noam Ross.
  • Explore rOpenSci projects
  • Do a deep dive into the rOpenSci blog history
  • Meet rOpenSci’s Executive Director, Noam Ross, and discuss how rOpenSci has evolved over the last 15 years.
• Tuesday July 7nd 2026, 09:00 Americas Pacific (16:00 UTC) “Debugging in R”, with Yanina Bellini Saibene and cohost Shannon Pileggi.
  • Read up on debugging in R
  • Meet community host, Shannon Pileggi, and discuss tips and tricks for debugging in R.

And remember, you can always cowork independently on work related to R, work on packages that tend to be neglected, or work on what ever you need to get done!

Software

New versions

The following ten packages have had an update since the last newsletter: cffr (v1.4.0), gutenbergr (v0.5.1), dataset (0.4.4), sofa (v0.4.1), saperlipopette (v1.0.0), medrxivr (snapshot), textreuse (v1.0.1), EDIutils (v2.1.0), rsi (v0.3.3), and UCSCXenaTools (v1.7.0).

The writexl package has a new maintainer, Bill Denney. NLMR is now maintained by Jakub Nowosad.

Software Peer Review

There are eighteen recently closed and active submissions and 4 submissions on hold. Issues are at different stages:

• One at ‘6/approved’:

  • saperlipopette, Create Example Git Messes. Submitted by Maëlle Salmon.

• Three at ‘5/awaiting-reviewer(s)-response’:

  • nycOpenData, Convenient Access to NYC Open Data API Endpoints. Submitted by Christian Martinez.

  • ernest, A Toolkit for Nested Sampling. Submitted by Kyle Dewsnap. (Stats).

  • pkgmatch, Find R Packages Matching Either Descriptions or Other R Packages. Submitted by mark padgham.

• Six at ‘4/review(s)-in-awaiting-changes’:

  • pvEBayes, Empirical Bayes Methods for Pharmacovigilance. Submitted by Yihao Tan. (Stats).

  • RAMEN, RAMEN: Regional Association of Methylome variability with the Exposome and geNome. Submitted by Erick Navarro-Delgado.

  • logolink, An Interface for Running NetLogo Simulations. Submitted by Daniel Vartanian.

  • rcrisp, Automate the Delineation of Urban River Spaces. Submitted by Claudiu Forgaci. (Stats).

  • priorsense, Prior Diagnostics and Sensitivity Analysis. Submitted by Noa Kallioinen. (Stats).

  • galamm, Generalized Additive Latent and Mixed Models. Submitted by Øystein Sørensen. (Stats).

• One at ‘3/reviewer(s)-assigned’:

  • ciecl, International Classification of Diseases ICD-10/ICD-11 for Chile. Submitted by Rodolfo Tasso.

• Five at ‘2/seeking-reviewer(s)’:

  • EpiStrainDynamics, Infer temporal trends of multiple pathogens. Submitted by Saras Windecker. (Stats).

  • lakefetch, Calculate Fetch and Wave Exposure for Lake Sampling Points. Submitted by jeremylfarrell.

  • RAQSAPI, A Simple Interface to the US EPA Air Quality System Data Mart API. Submitted by mccroweyclinton-EPA.

  • fcmconfr, Fuzzy Cognitive Map Analysis in R. Submitted by benroston. (Stats).

  • coevolve, Fit Bayesian Generalized Dynamic Phylogenetic Models using Stan. Submitted by Scott Claessens. (Stats).

• Two at ‘1/editor-checks’:

  • metasurvey, Reproducible Survey Data Processing with Step Pipelines. Submitted by Mauro Loprete.

  • LBDiscoverAnalysis, Co-occurrence Discovery Models and Visualization for Biomedical LBD. Submitted by Chao Liu.

Find out more about Software Peer Review and how to get involved.

On the blog

• New Mentoring Team with Experienced Mentors and New Voices by Alber Hamersson Sánchez Ipia, Pablo Paccioretti, Erick Navarro Delgado, Guadalupe Pascal, Andrea Gomez Vargas, Monika Avila Marquez, Luis D. Verde Arregoitia, Francisco Cardozo, Elio Campitelli, Pao Corrales, and Yanina Bellini Saibene. Meet the new team of mentors who will accompany our Champions in their projects for a year, sharing experience, guidance and enthusiasm. Other languages: Nuevo equipo de mentor(a|e)s con experiencia y nuevas voces (es).

• 15 Years of rOpenSci, and We’re Just Getting Started by The rOpenSci Team.

Calls for contributions

Calls for maintainers

If you’re interested in maintaining any of the R packages below, you might enjoy reading our blog post What Does It Mean to Maintain a Package?.

• charlatan, create fake data in R. Issue for volunteering.

• landscapetools, R package for some of the less-glamorous tasks involved in landscape analysis. Issue for volunteering.

• hddtools, Tools to discover hydrological data, accessing catalogues and databases from various data providers. Issue for volunteering.

Calls for contributions

Refer to our help wanted page – before opening a PR, we recommend asking in the issue whether help is still needed.

Package development corner

Some useful information for R package developers.

devtools’ install functions now officially deprecated

The devtools::install_ functions like devtools::install_github() have now been officially deprecated in favor of pak. The manual page contains a migration guide. Time to update older READMEs!

New Git commands!

New in Git: git history! Thanks to Hugo Gruson for bringing this to our attention.

Imagine…

Oops, that old commit’s message had a typo! git history reword

Oops, that old commit is too big! git history split (to split it into two commits)

Git release notes; Git docs.

roxygen2 8.0.0

A new version of roxygen2 hit CRAN. Among the highlights are improvements for R6, S7, inheritance, and better documentation of roxygen2 itself, including the vignette on extending roxygen2 with your own tags and roclets. Release announcement.

Comment DESCRIPTION files

Have you noticed this line in the release notes of R 4.6.0?

read.dcf() now recognizes lines starting with # as comment lines. By Dirk Eddelbuettel, Laurent Gatto and Hugo Gruson.

This means it’s becoming possible to add comments to DESCRIPTION files, for instance to note why you added a given dependency.

{cross} R package for running code with different package versions

Have you ever used Git worktree to load several R package versions at once (in different sessions), for instance to create a benchmark? Davis Vaughan made such endeavors much smoother with his cross R package. cross uses Git worktree under the hood to install different package versions in temporary libraries and then runs R code with these different versions. Furthermore, it has a user-friendly interface!

For example, to compare the released version of vctrs to the developement version on GitHub, you could run the following.

cross::bench_versions(pkgs = c("vctrs", "r-lib/vctrs"), {
 library(vctrs)
 x <- c(1, NA, 2, 3, NA)
 bench::mark(missing = vec_detect_missing(x))
})

On attribution with a footnote

James Balamuta wrote two thoughtful posts on a tricky situation and its resolution. TL;DR James’ work on portable R was used in another repository with only minimal acknowledgement through a footnote in a comment. However, those involved came together to work out a resolution through meaningful conversation.

Last words

Thanks for reading! If you want to get involved with rOpenSci, check out our Contributing Guide. This guide will help direct you to the right place, whether you want to make code contributions, non-code contributions, or contribute in other ways such as through sharing use cases. You can also support our work through donations.

If you haven’t subscribed to our newsletter yet, you can do so though our signup form. Until it’s time for our next newsletter, you can keep in touch with us through our website, Mastodon, or LinkedIn. See you soon!

To leave a comment for the author, please follow the link and comment on their blog: rOpenSci - open tools for open science.

---

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: rOpenSci News Digest, May 2026]]> 401509 https://www.r-bloggers.com/2026/05/survivor-now-includes-us50-and-au12/ Wed, 27 May 2026 09:54:21 +0000 http://gradientdescending.com/?p=3531 After a couple of big seasons survivoR v2.3.12 has been updated with US50 and AU12. Get it on: As of […]
The post survivoR now includes US50 and AU12 appeared first on Dan Oehm | Gradient Descending.

Continue reading: survivoR now includes US50 and AU12]]> [social4i size="small" align="align-left"]

After a couple of big seasons survivoR v2.3.12 has been updated with US50 and AU12.

Get it on:

• Github
• XLSX download
• JSON endpoints

As of writing, the survivoR package is closing in on 30,000 downloads, which is wild to me. Not only that, but the CRAN installs are probably the lowest volume of all formats. Installs directly from Git and the XLSX downloads are probably much higher. I never thought this data would be so popular.

I don’t have much of a social media presence, pretty much only Bluesky these days, and a very small following there. Regardless many people reach out to say thanks for pulling it together. It has enabled them to build their own websites (of which there are quite a few – Survivor-reference, survivorELO), be used for research papers, Master theses, school projects, help potential castaways with their Survivor auditions, help future Survivor players prep for their season, even then Survivordle web app.

I enjoy hearing about how others are utilizing it. If you have experience with it as well, please reach out—I am always keen to discuss Survivor!

The post survivoR now includes US50 and AU12 appeared first on Dan Oehm | Gradient Descending.

Introduction

In my previous job my work computer was a Windows desktop – yes, those were the days before laptops and hotdesking!

My PhD student was interested in Bayesian methods and we put together an R package which included some Stan models. I was always frustrated by how slowly these compiled on our Windows machines. A few years later, when I got a MacBook Air I was shocked how much faster they compiled.

On my Windows machine our mrbayes package takes 3 minutes 55 seconds to compile and install. On my M4 MacBook Air it takes 1 minute 16 seconds.

The following tips show how to improve those timings.

To generate the timings I used

time R CMD INSTALL --preclean .

Big win 1: Enable parallel compilations with the MAKEFLAGS environment variable

Set the MAKEFLAGS environment variable in your ~/.Renviron file. This controls how many make jobs run concurrently. Choose a number no larger than the number of processing cores your machine has. To find this run

# Windows - in a Git Bash shell
echo $NUMBER_OF_PROCESSORS

# macOS
sysctl -n hw.logicalcpu

# Ubuntu Linux
nproc

A reasonable starting point is your core count, or a few fewer to leave headroom for whatever else you’re doing during a compilation. For example,

# In ~/.Renviron
MAKEFLAGS=-j6

Close and restart R/RStudio after making this change.

On my Windows machine this reduced the build from 3:55 to 1:15. To find your own sweet spot empirically, see the example at the end of Big win 2.

Big win 2: Enable C/C++ compiler cache using ccache

Install ccache, I find it easiest to use a package manager, e.g.,

# macOS
brew install ccache

# Ubuntu/Debian Linux
apt install ccache

# Windows
winget install ccache

Whichever installation method you use make sure ccache is on your PATH after installation. You can test with, say,

ccache --version

To enable ccache, on macOS and Linux this goes in ~/.R/Makevars; on Windows it’s ~/.R/Makevars.win (create the directory and file if they don’t exist), set

# macOS
CC = ccache clang
CXX = ccache clang++
CXX17 = ccache clang++

# Windows and Linux
# Most Linux users will be on gcc by default
# Change to clang if you're using that
CC = ccache gcc
CXX = ccache g++
CXX17 = ccache g++

After a first compilation run for the cache to be generated, subsequent compilations are much faster.

• Windows, second compilation: 18 seconds
• M4 MacBook Air, second compilation: 5 seconds

Perhaps more importantly, if, say, your package has 5 models and you only amend the code for one of them, ccache knows to use the cache for the 4 unchanged models.

• Windows, second compilation, only 1 model edited: 1 minute 10 seconds
• M4 MacBook Air, second compilation, only 1 model edited: 19 seconds

You can verify ccache is working, by observing the timing decrease and by checking the output of

ccache -s

It is also useful to zero the ccache statistics before a timing run with

ccache -z

Testing which of your models takes the longest to compile

Here’s a quick script to test which model takes the longest to compile. Save it as say test.sh at the top level of your repo and add ^test\.sh$ to your .Rbuildignore file (to avoid an R CMD check NOTE about unknown files at the top level).

for model in inst/stan/*.stan; do
  cp "$model" "$model.bak"
  # Insert at the top of the file
  sed -i "1i // benchmark $(date +%s%N)" "$model"
  ccache -z
  SECONDS=0
  R CMD INSTALL --preclean . >/dev/null 2>&1
  echo "$(basename $model): ${SECONDS}s"
  ccache -s | grep -E "Hits|Misses" | head -2
  mv "$model.bak" "$model"
done

Finding your MAKEFLAGS sweet spot

With ccache installed you can now benchmark different -jN values cleanly (the ccache -C calls ensure each run is a cold compile, so you measure raw compilation cost rather than cache hits). You can increase the number sequence up to the number of processing cores your machine has.

for j in 1 2 3 4 6 8 10; do
  ccache -C >/dev/null
  echo "=== -j$j ==="
  SECONDS=0
  MAKEFLAGS=-j$j R CMD INSTALL --preclean . >/dev/null 2>&1
  echo "elapsed: ${SECONDS}s"
done

The timings on my MacBook Air were

=== -j1 ===
elapsed: 76s
=== -j2 ===
elapsed: 48s
=== -j3 ===
elapsed: 35s
=== -j4 ===
elapsed: 36s
=== -j6 ===
elapsed: 27s
=== -j8 ===
elapsed: 27s
=== -j10 ===
elapsed: 28s

My MacBook Air has 10 cores, but only 4 of those are performance cores, so I settled on -j6 as that is where my timings plateaued — and it leaves headroom for me inevitably checking my email during a compilation.

Big win 3: Combining these in GitHub Actions workflows

In my .github/workflows/R-CMD-check.yaml I have steps for these speedups. Firstly, to set MAKEFLAGS.

      - name: Set parallel compilation flags (Linux and macOS)
        if: runner.os != 'Windows'
        shell: bash
        run: |
          NCPUS=$(nproc 2>/dev/null || sysctl -n hw.logicalcpu)
          echo "Detected ${NCPUS} processors"
          echo "MAKEFLAGS=-j${NCPUS}" >> ~/.Renviron

      - name: Set parallel compilation flags (Windows)
        if: runner.os == 'Windows'
        shell: pwsh
        run: |
          Write-Output "Detected $env:NUMBER_OF_PROCESSORS processors"
          Add-Content -Path "$HOME\.Renviron" -Value "MAKEFLAGS=-j$env:NUMBER_OF_PROCESSORS"

You can also use ccache in GitHub Actions, as follows:

      # ccache speeds up Stan model compilation dramatically on warm cache.
      # Note: Windows support via ccache-action is documented as "probably works"
      # rather than fully stable; if it causes issues, scope this step to non-Windows.
      - name: Setup ccache
        uses: hendrikmuhs/ccache-action@v1.2.23
        with:
          # Key invalidates when Stan models or DESCRIPTION change.
          # Older caches partially seed new ones via restore-keys.
          key: ccache-${{ matrix.config.os }}-R-${{ matrix.config.r }}-${{ hashFiles('inst/stan/**/*.stan', 'DESCRIPTION') }}
          restore-keys: |
            ccache-${{ matrix.config.os }}-R-${{ matrix.config.r }}-
            ccache-${{ matrix.config.os }}-R-
          max-size: "2G"

      - name: Configure R to use ccache (Linux and macOS)
        if: runner.os != 'Windows'
        shell: bash
        run: |
          mkdir -p ~/.R
          if [ "$RUNNER_OS" = "macOS" ]; then
            cat >> ~/.R/Makevars <<'EOF'
          CC = ccache clang
          CXX = ccache clang++
          CXX14 = ccache clang++
          CXX17 = ccache clang++
          CXX20 = ccache clang++
          EOF
          else
            cat >> ~/.R/Makevars <<'EOF'
          CC = ccache gcc
          CXX = ccache g++
          CXX14 = ccache g++
          CXX17 = ccache g++
          CXX20 = ccache g++
          EOF
          fi
          echo "--- ~/.R/Makevars ---"
          cat ~/.R/Makevars

      - name: Configure R to use ccache (Windows)
        if: runner.os == 'Windows'
        shell: pwsh
        run: |
          New-Item -ItemType Directory -Force -Path "$HOME\.R" | Out-Null
          $makevars = @"
          CC = ccache gcc
          CXX = ccache g++
          CXX14 = ccache g++
          CXX17 = ccache g++
          CXX20 = ccache g++
          "@
          Add-Content -Path "$HOME\.R\Makevars.win" -Value $makevars
          Write-Output "--- ~/.R/Makevars.win ---"
          Get-Content "$HOME\.R\Makevars.win"

You can see the full file in my repo.

This reduced my ubuntu-latest run for r-release from 7 minutes 30 seconds to 4 minutes 49 seconds.

Big win 4: Switch to clang

I found that switching from gcc to clang gives a noticeable speedup; the single core compile time dropped from 3 minutes 55 seconds to 3 minutes flat on my Windows machine.

To do this you need to install clang. On Windows you install clang within RTools45 — more involved than on Linux, but doable.

# Windows within RTools45 Bash shell
# Launch C:\rtools45\ucrt64.exe
# You may need to close and reopen the shell after the first command
pacman -Syu
pacman -S mingw-w64-ucrt-x86_64-clang

# Ubuntu/Debian Linux
sudo apt install clang

At this point on Windows running

which clang

should return /ucrt/bin/clang.

Switch to clang in ~/.R/Makevars (if you’re not using ccache delete that prefix)

# On Linux
CC = ccache clang
CXX = ccache clang++
CXX14 = ccache clang++
CXX17 = ccache clang++
CXX20 = ccache clang++

and in ~/.R/Makevars.win on Windows

# On Windows
CC = ccache C:/rtools45/ucrt64/bin/clang.exe
CXX = ccache C:/rtools45/ucrt64/bin/clang++.exe
CXX17 = ccache C:/rtools45/ucrt64/bin/clang++.exe

Windows users will need to add the following to PATH

C:\rtools45\ucrt64\bin
C:\rtools45\usr\bin

You can verify things are working by running

R CMD config CXX17

I believe you need clang version 18 or later to see the speedups.

Small win 1: WSL users should use the native file system

Within WSL it is possible to access files from within its native Linux filesystem, i.e., within /home/user/..., and also on the Windows filesystem, e.g., in /mnt/c/.... I believe file operations are noticeably faster within /home/user/....

Naive guesses that made no difference

I had wondered whether running a non-debug compilation with say

pkgbuild::compile_dll(debug = FALSE)

would speed things up. It turns out it does not. For Stan models, most of the time is spent in C++ template instantiation by the compiler, not in optimisation passes — so disabling debug flags or lowering the optimisation level barely helps.

I also wondered whether using R on Windows Subsystem for Linux would speed things up just by virtue of being on Linux. It did not, timings using gcc on Windows and WSL Ubuntu were essentially identical. The advantage of using WSL is that it is easier to switch to using clang on Linux.

(Money no object) Big win 5: Switch to an Apple Silicon Mac

Apple silicon Macs have excellent single threaded performance, their unified memory architecture has very high bandwidth, they have large L1 and L2 caches, and fast NVMe SSDs. Together these produce very fast Stan model compilation times, even on the lowest end Apple Silicon Macs.

Summary

In summary, five big wins and one small win for speeding up Stan model compilation in R packages.

A specification tells a story.

Given describes the world before the action. When describes the action itself. Then describes what changed as a result.

If Given is the setup and Then is the payoff, When is where the plot happens, and how you write it determines whether your specifications stay readable for months or decay after the very first change you make to your project.

This article is part of a series on Behavior-Driven Development for Shiny applications. We’ve been building a data submission form from scratch, adding an email notification feature, and managing preconditions with Given steps.

Read the previous articles to get up to speed, or continue here to focus on how to write When steps.

1. Behavior-Driven Development in R Shiny: A Step-By-Step Example
2. Behavior-Driven Development in R Shiny: Setting Up Test Preconditions with Given Steps

Level-up your testing game! Grab your copy of the R testing roadmap.

The Purpose of When

When steps answer one question: What is the user doing?

Not what button they clicked. Not what input ID was set. What they did, from their perspective.

The difference matters because UI implementations change constantly, user behaviors change only when business requirements change. Input IDs get renamed, buttons get moved, forms get refactored into modules. But users still interact with the same form, they submit the same data. If your When steps name those interface details, they break every time something shifts. If they describe user intent, they survive almost any UI change.

A When step should read like a sentence from a user story.

Compare these two:

# Implementation-coupled! Breaks when the submit button is removed or chaned
when_i_click_submit_button <- function(context) {
  ...
}

# Behavior-focused! Survives any UI refactor
when_i_submit_entry_with_all_required_fields <- function(context) {
  ...
}

The second version tells you what the user accomplished, not how the app was wired.

That’s the goal.

Starting the App

When steps are where the Shiny app launches — but not directly.

As we established in the Given article, the app starts lazily: the first When step calls .start_driver(), which creates the driver with all configuration accumulated from Given steps. Every subsequent When step interacts with the already-running app.

#' tests/testthat/setup-dsl.R
when_i_submit_entry_with_all_required_fields <- function(context) {
  .start_driver(context)
  context$driver$fill_required_fields(
    url = "https://example.com/",
    title = "Example"
  )
  context
}

Calling .start_driver() in every When step is safe — it’s a no-op if the driver already exists:

.start_driver <- function(context) {
  if (is.null(context$driver)) {
    context$driver <- MyAppDriver$new(
      storage = context$storage,
      email_service = context$email_service,
      user = context$user
    )
  }
  context
}

This means Given steps are purely configuration — they never touch the running app. The app only starts when the user’s first action is about to happen. This makes Given steps fast, composable, and easy to reason about.

Implementing When Steps

The When step function itself should be thin. Its only job is to delegate to the driver and pass the context forward.

All implementation details live in the driver:

#' tests/testthat/setup-driver.R
MyAppDriver <- R6::R6Class(
  classname = "MyAppDriver",
  inherit = shinytest2::AppDriver,
  public = list(
    fill_required_fields = function(url, title) {
      self$set_inputs(url = url, title = title)
      self$click(input = "submit")
      self$wait_for_idle()
    },
    fill_all_fields = function(url, title, description) {
      self$set_inputs(url = url, title = title, description = description)
      self$click(input = "submit")
      self$wait_for_idle()
    }
  )
)

The driver encapsulates everything: which inputs to set, what to click, what to wait for. The When step just calls the right driver method by name.

This separation pays off immediately. If the submit button ID changes from "submit" to "btn_submit", you fix it in one place — the driver method — and every specification that calls fill_required_fields continues to work unchanged.

Note that the same driver method can be called from multiple When steps.

The driver is a translation layer. It translates user intent (fill required fields) into UI mechanics (set these inputs, click this button).

Naming When Steps

Good names follow a simple pattern: when_<who>_<verb>_<object>.

The verb describes the action from the user’s point of view. The object is what they acted on — expressed in domain language, not UI language.

# Good! We use domain language
when_i_submit_entry_with_all_required_fields()
when_i_submit_entry_with_missing_title()
when_i_delete_the_first_entry()
when_i_inspect_the_submissions_list()

# Bad! UI details are leaking in
when_i_click_the_submit_button()
when_i_clear_the_title_input()
when_i_navigate_to_slash_list()

Domain names survive redesigns. UI names couple your specifications to the current implementation.

One more rule: one action per When step. A step named when_i_fill_and_submit_the_form() is doing two things. That makes it harder to compose and harder to understand which action caused which outcome.

# One action per step — composable
when_i_fill_in_the_form(url = "https://example.com", title = "Example") |>
  when_i_submit_the_form()

# Doing too much — hard to reuse
when_i_fill_and_submit_the_form()

Smaller, focused steps let you test multi-step flows precisely.

My Thoughts on Navigation Steps

You might think, this is an easy example.

Real apps aren’t single-page forms. They have navigation, multiple pages, modals, they might feel like navigating through a maze. It’s not that simple to interact with all the features! It’s not, but we can apply the same principles to more complex interactions.

Here’s what I do.

In apps we have pages. Users need to navigate to those pages before they can interact with them. How do we name When steps that represent navigation?

I prefer to name them based on their intent, what they’re trying to achieve, rather than the mechanics of how they get there. For example:

# Good! Focuses on the user's goal, not the UI structure
when_i_inspect_submissions()
when_i_open_the_first_entry()
when_i_change_settings(...) # composite step, implicitly navigates to the settings page and changes settings

# Bad! Reveals that there is a "submissions" page, what if we move submissions listing to another page?
when_i_navigate_to_submissions_page()

This way, the step names reflect the user’s goals, not the UI structure. If we later change how navigation works — say we add another page and move our targeted component there — we can update the driver methods without changing the When step names.

That’s how specifications stay stable and focused on user behavior, not implementation details.

That’s how specifications stay stable and don’t need updating every time the UI changes.

Multi-Step Flows

When a scenario involves several user actions in sequence, chain them:

it("should allow editing a submitted entry", {
  given_existing_entries() |>
    given_an_authenticated_user() |>
    when_i_open_the_first_entry() |>
    when_i_update_the_title(title = "Updated Title") |>
    when_i_submit_the_form() |>
    then_the_entry_has_title("Updated Title")
})

Each step does one thing. The sequence reads like a user story. And if the assertion fails, you can see exactly which action preceded the broken state.

This is also where the context pipeline pattern earns its keep. Every step receives context and returns context, so you can add or remove steps in the middle without restructuring the whole test:

it("should show confirmation after submitting", {
  given_no_content() |>
    given_an_authenticated_user() |>
    when_i_submit_entry_with_all_required_fields() |>
    then_entry_is_submitted()
})

The pipeline works because every step function — Given, When, or Then — has the same contract: receive context, return context. That uniformity is what lets pipe operator connect them.

Here is the pattern on its own, with no Shiny involved. The When steps are thin — each one just calls a driver method and returns context:

# A minimal driver that records what it was asked to do
MockDriver <- function() {
  form <- NULL
  submitted <- FALSE
  list(
    fill_required_fields = function(url, title) {
      form <<- list(url = url, title = title)
    },
    submit = function() {
      submitted <<- TRUE
    },
    get_form = function() form,
    was_submitted = function() submitted
  )
}

# When steps are thin: call a driver method, return context
when_i_fill_in_the_form <- function(context, url, title) {
  context$driver$fill_required_fields(url = url, title = title)
  context
}

when_i_submit_the_form <- function(context) {
  context$driver$submit()
  context
}

ctx <- new.env(parent = emptyenv())
ctx$driver <- MockDriver()

ctx |>
  when_i_fill_in_the_form(url = "https://example.com", title = "Example") |>
  when_i_submit_the_form()

The step functions contain no UI details — those live entirely inside MockDriver. Swapping the driver is all it takes to point the same steps at a different implementation.

Testing Error Paths

User interactions that trigger failures are still described with When steps. The step describes what the user attempted, not whether it succeeded.

it("should require all required fields", {
  given_no_content() |>
    given_an_authenticated_user() |>
    when_i_submit_entry_with_missing_required_fields() |>
    then_i_am_prompted_to_provide_required_fields()
})

The When step looks identical to the happy-path version — which is the point. The user tried to submit with missing fields. The difference lives in the Then step, which checks what the app communicated back to the user.

when_i_submit_entry_with_missing_required_fields <- function(context) {
  .start_driver(context)
  context$driver$fill_required_fields(
    url = "https://example.com/",
    title = ""  # Missing required field
  )
  context
}

You don’t need to create special “error path” When steps with names like when_i_incorrectly_submit(). Name them after what the user did, not after whether it was correct. The Then step will check whether the app responded correctly to that action.

Parameterising When Steps

Sometimes you want to test multiple variations of the same action. Rather than creating a separate step for each variation, make the step accept parameters:

when_i_submit_entry_with_title <- function(context, title) {
  .start_driver(context)
  context$driver$fill_required_fields(
    url = "https://example.com/",
    title = title
  )
  context
}

Now you can express different scenarios clearly:

it("should reject titles that are too short", {
  given_no_content() |>
    given_an_authenticated_user() |>
    when_i_submit_entry_with_title(title = "ab") |>
    then_i_am_prompted_to_provide_required_fields()
})

it("should accept titles of the minimum length", {
  given_no_content() |>
    given_an_authenticated_user() |>
    when_i_submit_entry_with_title(title = "abc") |>
    then_there_are_entries(n = 1)
})

Use parameters when the difference between scenarios is data, not behavior. Keep the step name stable; let the arguments express the variation.

The Cucumber Way

cucumber allows you to write multiple variants of scenarios with Scenario Outline and Examples keywords.

In a cucumber feature file this scenario would look like this:

Feature: Data submission form validation
  Scenario Outline: Submitting entries with different title lengths
    Given no content
    And an authenticated user
    When I submit an entry with title <title>
    Then I am prompted to provide required fields
    Examples:
      | title |
      | "ab"  |
      | "abc" |

In our implementation we don’t have a built-in way to express scenario outlines, but we can achieve the same effect by writing multiple scenarios that call the same parameterised When step.

Signs a When Step Is Too Complex

A When step that runs more than a few lines of logic is usually trying to do too much.

1. If a step is setting state — checking whether something exists before interacting with it, or conditionally taking different actions — that logic belongs in Given steps or in the driver implementation, not the When step itself.
2. If a step has branching logic — doing one thing when a flag is set and another when it isn’t — split it into two steps with clear names.
3. If a step is making assertions — checking that something happened mid-flow — move those checks to Then steps. A When step should interact and nothing more. Mixing assertions into When steps makes failures ambiguous: you can’t tell whether the setup was wrong or the outcome was wrong. There might be exceptions to assert whether we’re on the right path in a multi-step flow – to fail eagerly – but generally, try to keep assertions in Then steps.

The Pattern in Practice

Here’s the full shape of a specification using the When steps we’ve built:

describe("data submission", {
  it("should submit entry and send notification", {
    given_no_content() |>
      given_an_authenticated_user() |>
      given_email_service_is_available() |>
      when_i_submit_entry_with_all_required_fields() |>
      then_entry_is_submitted() |>
      then_email_notification_is_sent()
  })

  it("should require all required fields", {
    given_no_content() |>
      given_an_authenticated_user() |>
      when_i_submit_entry_with_missing_required_fields() |>
      then_i_am_prompted_to_provide_required_fields()
  })

  it("should allow editing a submitted entry", {
    given_existing_entries() |>
      given_an_authenticated_user() |>
      when_i_open_the_first_entry() |>
      when_i_update_the_title(title = "Updated Title") |>
      when_i_submit_the_form() |>
      then_the_entry_has_title("Updated Title")
  })
})

The specifications read naturally. None of them mention input IDs, button selectors, or shinytest2. Those details live in the driver, hidden away from the scenarios they serve.

Note how specifications don’t even reveal it’s a Shiny app. They don’t reveal how we interact with it — shinytest2, selenider, anything goes. It’s all hidden and irrelevant to the story we’re telling. You can swap the internals and keep the same specifications true. You can even swap the system implementation entirely, and as long as it still satisfies the same business requirements, the specifications don’t need touching. It’s live and executable — not a list of tickets in an external system no one reads.

That’s exactly the opposite of what happens when we write tests that are tightly coupled to the UI structure.

Wrapping up

When steps model user behavior, not UI implementation.

That distinction is what keeps specifications readable days, weeks, months later, when the form has been refactored, button IDs have changed, and nobody remembers what the original implementation looked like. Keep steps focused on one action. Put the implementation details in the driver. Use domain language, not UI language.

Specifications written this way tend to survive refactors you haven’t thought of yet.

Reposted from the original at https://blog.stephenturner.us/p/ctrlvee-extract-external-r-code-insert-inline-positron-rstudio-addin. 

Ever find yourself looking through a pkgdown page or a Quarto book, copying and pasting code chunks from your browser into your IDE? I do, and it’s a minor annoyance.1

My friend and colleague VP Nagraj published a new R package called ctrlvee that makes this a lot easier.

• CRAN: https://cran.r-project.org/package=ctrlvee

• GitHub: https://github.com/vpnagraj/ctrlvee

It does one thing. Put your cursor anywhere in an R script in Positron or RStudio, call the add-in, provide a URL, and a few milliseconds later you’ll have all the code from that page in your editor, separated by chunk boundaries (along with some metadata and a note to check the license!).

The package README provides a demonstration using the “Data Validation and QA” chapter of my Data Science Team Training book (dstt.stephenturner.us).

1. Install the package: install.packages("ctrlvee")

2. Run the add-in. In Positron you’ll open the command palette, search for Run RStudio Addin, then extract external R code and insert inline. You’ll get a modal asking you for a URL.

3. Paste one in. E.g., https://dstt.stephenturner.us/validation.html

4. The R code from the website appears in your editor

Here’s a demo.

Here’s what the extracted/inserted code looks like, from this source.

# -----------------------------------------------------------------
# Chunks fetched by ctrlvee from: https://dstt.stephenturner.us/validation.html
# Strategy: Rendered HTML page
# Date: 2026-05-16 05:14:44
# Chunks: 8
# NOTE: Check the source license before reusing this code.
# -----------------------------------------------------------------

flu <- data.frame(
    week = c(1, 2, 3, 4, 4),
    county = c("Fairfax", "Arlington", NA, "Loudoun", "Loudoun"),
    disease = c("Flu", "Flu", "Flu", "Flu", "Flu"),
    cases = c(23, 41, 18, -5, 12),
    rate = c(2.1, 3.8, 1.6, NA, 1.1)
)

flu

# ---- chunk boundary ----

if (any(flu$cases < 0, na.rm = TRUE)) {
    stop("Negative case counts detected. Inspect raw data before proceeding.")
}

# ---- chunk boundary ----

stopifnot(
    "Negative case counts" = all(flu$cases >= 0, na.rm = TRUE),
    "Missing county values" = !anyNA(flu$county),
    "Duplicate records" = !anyDuplicated(flu[, c("week", "county")])
)

# ---- chunk boundary ----

install.packages("pointblank")

# ---- chunk boundary ----

library(pointblank)

agent <- create_agent(tbl = flu, label = "Weekly flu surveillance") |>
    col_vals_gte(
        columns = cases,
        value = 0,
        label = "Case counts must be non-negative"
    ) |>
    col_vals_not_null(
        columns = c(week, county),
        label = "Week and county cannot be missing"
    ) |>
    rows_distinct(
        columns = c(week, county),
        label = "No duplicate week/county records"
    ) |>
    interrogate()

agent

# ---- chunk boundary ----

create_agent(tbl = flu, label = "Weekly flu surveillance — extended") |>
    col_is_numeric(
        columns = c(cases, rate),
        label = "Case count and rate must be numeric"
    ) |>
    col_vals_in_set(
        columns = disease,
        set = c("Flu", "COVID-19", "RSV"),
        label = "Disease must be from the approved list"
    ) |>
    col_vals_between(
        columns = week,
        left = 1,
        right = 52,
        label = "Week must be between 1 and 52"
    ) |>
    col_vals_gte(
        columns = rate,
        value = 0,
        na_pass = TRUE,
        label = "Rate must be non-negative (NAs allowed)"
    ) |>
    interrogate()

# ---- chunk boundary ----

if (!all_passed(agent)) {
    stop("Data validation failed. Review the agent report before proceeding.")
}

# ---- chunk boundary ----

library(readr)
library(pointblank)

flu <- read_csv("data/flu-2024.csv")

# Validate immediately after reading
agent <- create_agent(tbl = flu, label = "flu-2024 validation") |>
    col_vals_gte(columns = cases, value = 0, label = "No negative counts") |>
    col_vals_not_null(columns = c(week, county), label = "No missing keys") |>
    rows_distinct(columns = c(week, county), label = "No duplicate records") |>
    interrogate()

if (!all_passed(agent)) {
    stop("Validation failed — see agent report above.")
}

Adding patterns to plots is a great way to improve accessibility (making plots colorblind-friendly) and to add an extra dimension of information. The ggpattern package provides a rich set of tools to achieve this in ggplot2.

# Prepare data
df <- mtcars %>%
  mutate(
    cyl = factor(cyl),
    gear = factor(gear),
    am = factor(am, labels = c("Automatic", "Manual"))
  )

# Create barplot
ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    color = "black",
    pattern_fill = "black",
    pattern_spacing = 0.02,
    pattern_key_scale_factor = 1.5
  ) +
  scale_pattern_manual(values = c("stripe", "crosshatch")) +
  labs(
    title = "MPG by Cylinders, Gear, and Transmission",
    subtitle = "Patterns representing transmission type",
    x = "Cylinders",
    y = "Miles Per Gallon (MPG)",
    fill = "Gears",
    pattern = "Transmission"
  ) +
  theme_minimal()

p_base <- ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    color = "black",
    pattern_spacing = 0.02
  ) +
  theme_minimal() +
  labs(title = "Default vs Scaled Legend")

# Default legend
p1 <- p_base + labs(subtitle = "Default (key_scale_factor = 1)")

# Scaled legend
p2 <- p_base + 
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    color = "black",
    pattern_spacing = 0.02,
    pattern_key_scale_factor = 0.5
  ) +
  labs(subtitle = "Scaled (key_scale_factor = 0.5)")

p1

p2

# Spacing effect
ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    pattern_spacing = 0.2, # Wider spacing
    pattern_key_scale_factor = 1.5
  ) +
  labs(title = "Effect of pattern_spacing = 0.5 (Sparser)") +
  theme_minimal()

# Density effect
ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    pattern_density = 0.8, # Higher density
    pattern_key_scale_factor = 1.5
  ) +
  labs(title = "Effect of pattern_density = 0.8 (Denser)") +
  theme_minimal()

ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern = am),
    stat = "identity",
    position = "dodge",
    pattern_fill = "white", # White pattern lines
    pattern_density = 0.5, # make it bigger so that the color is easier to see
    pattern_spacing = 0.02,
    pattern_key_scale_factor = 1.5
  ) +
  labs(title = "Effect of pattern_fill = 'white'") +
  theme_minimal()

ggplot(df, aes(x = cyl, y = mpg, fill = gear)) +
  geom_bar_pattern(
    aes(pattern_shape=gear),
    pattern = "regular_polygon",
    stat = "identity",
    position = "dodge",
    pattern_size = 0.5, # stroke line width
    pattern_density = 0.5, # increase the shape size
    pattern_spacing = 0.03,
    pattern_key_scale_factor = 1.5
  ) +
  scale_pattern_shape_manual(values = c("3"="circle", "4"="convex6", "5"="square")) + # manually set shapes
  labs(title = "Effect of mapping pattern_shape to gear") +
  theme_minimal()

sudo apt-get install libudunits2-dev
sudo apt install libproj22
sudo apt-get install libgdal-dev

If you’ve read any of my past posts you know I like to program in several different languages, some of which I like more than others. Sometimes a problem calls for a particular language to be used, and with that comes adjusting one’s brain to thinking in that language and using the appropriate idioms to leverage that language’s features. But what if I don’t want to?

I don’t want to

The line between R and Python has been heavily blurred the last few years, particularly with {reticulate} enabling us to use Python within R code, RStudio rebranding as Posit and taking on a strong Python development effort, releasing Positron as a multi-language IDE, and Quarto being a multi-language rethink of Rmarkdown.

I occasionally need to use Python directly – an SDK wrapping an API exists and I don’t particularly want to spend a lot of time writing my own R version, especially before I know what I want to get out of the endpoints. At this point I tend to bump up against my muscle-memory from R and try to use functions I’m familiar with from R, but which don’t actually exist in Python. Now, that might sometimes be because the pattern I’m trying to encode simply has a different name in Python; instead of an sapply(x, f)

sapply(c(2, 3, 4, 5), \(x) x ^ 2)
## [1]  4  9 16 25

I should reach for map, in which case I am reminded that this produces a lazy iterator that doesn’t show me the results

map(lambda x: x ** 2, [2, 3, 4, 5])
## <map object at 0x10d7fbee0>

and so I need to wrap it into a list to get the values out

list(map(lambda x: x ** 2, [2, 3, 4, 5]))
## [4, 9, 16, 25]

Or, I could use a list comprehension which isn’t lazy

[v ** 2 for v in [2, 3, 4, 5]]
## [4, 9, 16, 25]

That’s the idiom that I should be reaching for. Sure.

Other times there’s a package I need to use and a slightly different way of approaching the problem. In R I love the table() function for getting histogram-like counts of the unique values of a vector

table(c("b", "a", "c", "a", "b", "a"))
## 
## a b c 
## 3 2 1

which in Python looks like

from collections import Counter

sorted(Counter(["b", "a", "c", "a", "b", "a"]).items())
## [('a', 3), ('b', 2), ('c', 1)]

Probably Pythonistas remember that idiom and the package to import and the .items() extractor and the fact that they maybe want to sort the result. But I kept coming back to a question I ask myself: what if I don’t want to? Why is there not a function that wraps this idiom? If there was, why not just call it “table”? Admittedly, it’s far from the catchiest, most memorable, or most useful name, but it’s immediately recognisable to an R user (ditto for “sapply”).

One approach I considered here was to just call R from Python. That can be done, but I doubt I or anyone else wants to deal with that every time we want to iterate over a list. There’s a package on the Python package index which seems to support this nicely: https://pypi.org/project/r-functions/ but it’s wrappers around individual R files, via RScript. I’m thinking more along the lines of ‘native Python with an R interface’.

Python is an object-oriented language, but it has functions, so why not make one

from collections import Counter

def table(x):
    return dict(sorted(Counter(x).items()))

table(["b", "a", "c", "a", "b", "a"])
## {'a': 3, 'b': 2, 'c': 1}
def sapply(x, func):
    return [func(v) for v in x]
  
sapply([2, 3, 4, 5], lambda x: x ** 2)
## [4, 9, 16, 25]

and have a nicer function interface to apply these idioms? I thought about this a bit longer, and realised there’s lots of functions I use in R that I wish I could use in Python. An idiom for finding the index of elements of a ‘vector’ (list in Python) which are true (TRUE in R, True in Python) is

[i for i, v in enumerate(x) if v]

but I just want to call which(x)

which(c(FALSE, FALSE, TRUE, FALSE , TRUE))
## [1] 3 5

so why not define this

def which(x):
    return [i for i, v in enumerate(x) if v]
  
which([False, False, True, False, True])
## [2, 4]

(remembering that Python is 0-indexed).

How far could one take this? Quite a long way!

I thought more about what differences would need to be accounted for, and one that immediately came to mind was that R is vectorised. If I was to recreate R’s character counting function nchar(s) as essentially len(s), I’d need to consider whether I wanted it to work on a single string or a ‘vector’ of strings

In R:

nchar(c("these", "all", "have", "different", "lengths"))
## [1] 5 3 4 9 7

But in Python, len() expects a single value, so it calculates the length of the list

len(["these", "all", "have", "different", "lengths"])
## 5

The ‘proper’ way to do it is to map over the list

[len(s) for s in ["these", "all", "have", "different", "lengths"]]
## [5, 3, 4, 9, 7]

but again, why do I need to use an idiom for this? What if I just made a decorator to change a regular function to a vectorised one by applying this list comprehension internally when it’s passed a list (or a tuple), and which otherwise just evaluates the function with the argument?

import functools

def make_vec(func):
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        if isinstance(args[0], (list, tuple)):
            return [func(xi, *args[1:], **kwargs) for xi in args[0]]
        return func(*args, **kwargs)
    return wrapper

@make_vec
def my_len(s):
    return len(s)

my_len(["these", "all", "have", "different", "lengths"])
## [5, 3, 4, 9, 7]

and I could name it… “nchar”!

The other use-case that came to mind was Elio venting (and referencing a post to which I also wrote a sort of response) that they needed to list the files in the current directory

Post by @eliocamp@mastodon.social

View on Mastodon

with the idiom

import os

[os.path.join(path, f) for f in os.listdir(path)]

The supplied suggestions included

from pathlib import Path

list(Path(path).iterdir())

(just rolls off the tongue, doesn’t it?) which returns a list of PosixPath() objects and is hardly easy to parse visually.

So, why not have a function?!?

import os

def list_files(path):
    return [os.path.join(path, f) for f in os.listdir(path)]

path = "path/to/files"

list_files(path)
## ['path/to/files/file1.txt', 'path/to/files/file2.txt', 'path/to/files/file3.csv']

I would have liked to call this list.files() but, since Python strictly uses the dot for method calling, it can’t be that.

This then raises the question of “should I support the arguments already in the R functions?” In this case, should it support a recursive argument? Yes, that adds complexity, but it’s surely do-able. At this point I reached for some AI assistance and had Claude help me to implement as many functions as we could think of, supporting as many common arguments as possible. This involved extending the decorator to support vectorising other arguments (which also need to be careful about dots).

On testing it out, it looked like we had something viable.

One last piece I wanted to support, though: the which() example above extracts the elements of a logical vector which are True, but in order to build that vector in the first place, I would naturally leverage R’s vectorisation as an array language. The two steps involved here are to first compute the comparison resulting in a logical vector, then to use which() to identify the indices of those which are true

which(c("c", "b", "a", "c", "a", "b") == "a")
## [1] 3 5

The vectorisation decorator above doesn’t help here, because it’s at the point of == that we want to vectorise

['c', 'b', 'a', 'c', 'a', 'b'] == 'a'
## False

This is False because the character 'a' is not equal to the given list.

The appropriate idiom is once again to use a list comprehension

which(x == 'a' for x in ['c', 'b', 'a', 'c', 'a', 'b'])
## [2, 4]

The solution I’m fond of is to create a new ‘Vec’ class which wraps binary operators with a list comprehension, again abstracting away this detail. This means implementing __eq__, __add__, __and__ and lots of other binary operations, but with that, and a wrapper to create such an object, the comparison operators can be vectorised

vals = vec(['c', 'b', 'a', 'c', 'a', 'b'])
which(vals == 'a')
## [2, 4]

Not pristine, but quite clean, if you ask me.

With all these pieces in place, adding implementations for common base R functions including most arguments and a way to vectorise lists, I wrapped everything up into a Python package (my first) to learn how to do it.

The workflow isn’t particularly painful, with my biggest complication being different versions of Python supporting different requirements in pyproject.toml, and so some GitHub Actions are failing because of that.

As part of building out the implementations I had Claude add tests for each of the functions with some expected values – if I do want to improve some of the idioms internally, I want to ensure I don’t change the values produced. That works for having any testing at all, but how can I be sure that I’m reproducing what I would get if I was working in R? One option was to just run all of the test functions by hand and confirm that the values look similar enough, accounting for list vs vector and 0 vs 1 indexing. Instead, Claude managed to write an adaptor for pytest which does the realignment of e.g. list_files to list.files (and similarly for arguments), realigns the indexing where needed, and runs all existing tests directly in R via rpy2 (skipping over some for which I don’t have tests yet). I’m disabling automated testing of this because I suspect it could get flaky dealing with both R and Python on GitHub Actions, but I can confirm that all the current tests pass.

I wanted to have a documentation website similar to what we have via {pkgdown} and came across quartodoc which is what the Python version of {pins} uses. Getting that to work required downgrading a specific Python dependency, but was otherwise painless.

I have a working package locally – how do I share it? This seemed like the perfect opportunity to learn what the release process looks like for Python. I have a handful of packages on CRAN and one on Bioconductor, and the process there is far from frictionless, with the side-effect that there’s some trust you can place on the interoperability of packages and minimal (automated) code checking. While Python is more ‘wild west’ in terms of what can be uploaded, it’s really nice to see that they do have an entirely separate test server where you can upload your package and see how it looks. I’m reminded of the quote

Everybody has a testing environment. Some people are lucky enough to have a totally separate environment to run production in.

Given that it’s not currently possible to run 100% of the CRAN checks locally (and even some that you can give a different result to what’s on their systems) this does make me a little jealous. I wonder whether the decrease in load from rejecting failing submissions would offset supporting a test submission server.

All went well pushing to the test server (via an authentication key) and I managed to build up the courage to push to the production instance… it’s live!

rfuns logo – R functions in Python… are fun

and the documentation site isn’t too bad, either (in my opinion).

This means that you can now run

uv add rfuns

(or the equivalent in whatever virtual environment management configuration you’re using, e.g. pip install rfuns) and start using some R functions directly in Python!

Depending on how you like to manage your imports, you can import everything

from rfuns import *

which([False, False, True, False, True])
## [2, 4]

or, if you prefer to namespace

import rfuns as r

r.which([False, False, True, False, True])
## [2, 4]

The list of functions currently imported, grouped into sections is:

nchar(["these", "all", "have", "different", "lengths"])
## [5, 3, 4, 9, 7]
grepl("ar", ["frog", "carpet", "basket", "dart"])
## [False, True, False, True]
sqrt([36, 81, 9])
## [6.0, 9.0, 3.0]

seq(5)
## [0, 1, 2, 3, 4]
seq(from_=0, to=10, by=2)
## [0, 2, 4, 6, 8, 10]

setdiff([5, 2, 4, 1], [2, 1])
## [5, 4]

set([5, 2, 4, 1]) -  set([2, 1])
## {4, 5}

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.5.3 (2026-03-11)
##  os       macOS Tahoe 26.3.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Australia/Adelaide
##  date     2026-05-22
##  pandoc   3.6.3 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
##  quarto   1.7.31 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date (UTC) lib source
##  blogdown      1.23    2026-01-18 [1] CRAN (R 4.5.2)
##  bookdown      0.46    2025-12-05 [1] CRAN (R 4.5.2)
##  bslib         0.10.0  2026-01-26 [1] CRAN (R 4.5.2)
##  cachem        1.1.0   2024-05-16 [1] CRAN (R 4.5.0)
##  cli           3.6.5   2025-04-23 [1] CRAN (R 4.5.0)
##  devtools      2.4.6   2025-10-03 [1] CRAN (R 4.5.0)
##  digest        0.6.39  2025-11-19 [1] CRAN (R 4.5.2)
##  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.5.0)
##  evaluate      1.0.5   2025-08-27 [1] CRAN (R 4.5.0)
##  fastmap       1.2.0   2024-05-15 [1] CRAN (R 4.5.0)
##  fs            1.6.7   2026-03-06 [1] CRAN (R 4.5.2)
##  glue          1.8.1   2026-04-17 [1] CRAN (R 4.5.2)
##  htmltools     0.5.9   2025-12-04 [1] CRAN (R 4.5.2)
##  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.5.0)
##  jsonlite      2.0.0   2025-03-27 [1] CRAN (R 4.5.0)
##  knitr         1.51    2025-12-20 [1] CRAN (R 4.5.2)
##  lattice       0.22-9  2026-02-09 [1] CRAN (R 4.5.3)
##  lifecycle     1.0.5   2026-01-08 [1] CRAN (R 4.5.2)
##  magrittr      2.0.4   2025-09-12 [1] CRAN (R 4.5.0)
##  Matrix        1.7-4   2025-08-28 [1] CRAN (R 4.5.3)
##  memoise       2.0.1   2021-11-26 [1] CRAN (R 4.5.0)
##  otel          0.2.0   2025-08-29 [1] CRAN (R 4.5.0)
##  pkgbuild      1.4.8   2025-05-26 [1] CRAN (R 4.5.0)
##  pkgload       1.5.0   2026-02-03 [1] CRAN (R 4.5.2)
##  png           0.1-9   2026-03-15 [1] CRAN (R 4.5.2)
##  purrr         1.2.2   2026-04-10 [1] CRAN (R 4.5.2)
##  R6            2.6.1   2025-02-15 [1] CRAN (R 4.5.0)
##  Rcpp          1.1.1   2026-01-10 [1] CRAN (R 4.5.2)
##  remotes       2.5.0   2024-03-17 [1] CRAN (R 4.5.0)
##  reticulate    1.45.0  2026-02-13 [1] CRAN (R 4.5.2)
##  rlang         1.1.7   2026-01-09 [1] CRAN (R 4.5.2)
##  rmarkdown     2.30    2025-09-28 [1] CRAN (R 4.5.0)
##  rstudioapi    0.18.0  2026-01-16 [1] CRAN (R 4.5.2)
##  sass          0.4.10  2025-04-11 [1] CRAN (R 4.5.0)
##  sessioninfo   1.2.3   2025-02-05 [1] CRAN (R 4.5.0)
##  usethis       3.2.1   2025-09-06 [1] CRAN (R 4.5.0)
##  vctrs         0.7.1   2026-01-23 [1] CRAN (R 4.5.2)
##  withr         3.0.2   2024-10-28 [1] CRAN (R 4.5.0)
##  xfun          0.56    2026-01-18 [1] CRAN (R 4.5.2)
##  yaml          2.3.12  2025-12-10 [1] CRAN (R 4.5.2)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## 
## ─ Python configuration ───────────────────────────────────────────────────────
##  python:         /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/bin/python
##  libpython:      /Users/jono/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/libpython3.12.dylib
##  pythonhome:     /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV:/Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV
##  virtualenv:     /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/bin/activate_this.py
##  version:        3.12.12 (main, Oct 28 2025, 11:52:25) [Clang 20.1.4 ]
##  numpy:          /Users/jono/.cache/uv/archive-v0/Q1veGTfRq3GBaNYBXjagV/lib/python3.12/site-packages/numpy
##  numpy_version:  2.4.6
##  
##  NOTE: Python version was forced by VIRTUAL_ENV
## 
## ──────────────────────────────────────────────────────────────────────────────

Over at Daring Fireball, John Gruber makes a passing observation about the Apple Sports app:

I’ve got some gripes about certain specific aspects of Apple Sports. Like, where does one even start to explain how much is wrong with their zero-sum visualization of team stats? Has anyone ever even seen a presentation like that before? Anyone?

That “Anyone” link lands over here. Hi everyone! The team stats image is quite confusing. It’s a summary of a game between the San Antonio Spurs and the Oklahoma City Thunder. I don’t know much about basketball, but I do know a bit about data visualization and in a pleasing coincidence my former student Josh Fink is the A-VP of Basketball Data Science for the Spurs. Here is the image that John objected to:

I just finished driving a very long way up the side of the country, so I’m kind of tired. But even allowing for that, boy, this way of representing things really is quite confusing. Not being an Apple Sports user I had to look at it for a bit to understand what was happening. But, now that it has given me a headache, I can kind of see why whoever designed this ended up in the undoubtedly bad place they did.

Before I get to why I have some sympathy for the designer, why did I find this representation of these numbers so disorienting? It’s not just just because I’ve been driving for nine hours. John is right to call the picture a “Zero Sum” representation. The design strongly suggests to the viewer that, within each row, we’re looking at each team’s share of a total. Each pair of black and blue lines seem to be vying for control of their whole row, with the longest line being the “winner” in each case.

This sort of representation would make perfect sense for a measure that really was zero sum. Take an example from a properly good sport, like rugby. There, like in basketball, to a first approximation a team either has the ball or it doesn’t.¹ But there’s no shot clock in rugby, and possession routinely gets turned over without the game stopping. So, knowing that Team A had 65% possession is not only informative, it also immediately entails that Team B had 35%. You could show that with a representation like one of the rows above.

Literally none of the measures in the Basketball data above are zero-sum in this way. Both teams could shoot 100% from the free throw line, or zero percent. But because the first three measures shown are percentages, this reinforces the zero-sum impression given by the lines. It certainly did that in my case. But then, starting with Assists, the remaining rows are just absolute numbers. When I started looking at the absolute numbers, I got confused a second time by the length of the lines. “Oh so it’s not a share, it’s the value” I thought—but no, they do correspond in terms of relative proportions to the teams share within each row. But they’re not really shares they’re just magnitudes. But they have to be shown in a fixed space and we want to make them relatively comparable somehow so … Argh.

It would be nice if there were One Weird Trick to fully fix this figure. But I’m not sure that there is. For example, at a minimum we could redraw these numbers to reflect the fact that they’re not zero-sum. Keep each measure as a row (i.e. on the y-axis) but have the lines, or columns, be side by side within each category instead of facing off. Like this:

This view at least lets you immediately see who “won” each measure. The viewer can just directly compare the length of the bars in each category. People are really good at doing that accurately. In that sense it’s much less confusing than the original. But there’s still a lot wrong with it. The core problem is that when we draw a graph like this, we’re usually putting the same kind of thing (e.g. countries, or religious groups, or sports teams) on the y-axis, and then seeing how different their scores are on some single measure (e.g. GDP, or number of adherents, or average points scored per game), which we put on the x-axis. Maybe we use color to break things out by some third measure as well.² In this case, I’ve just labeled the x-axis as generically as possible. “Value” covers the range of all the measures. The lowest value is 5, in Largest Lead. The highest is 88, in Free Throw %. But these numbers are not meaningfully comparable. The graph encourages us to compare across as well as within categories. But while within-category comparisons are meaningful, the between-category ones are not. There were way more Bench Points than Blocks in the game. But that is not a useful thing to know.

Knowing who won each measure isn’t nothing. It can be informative about how the game went, maybe especially when a team won the game but “lost” on a number of the measures. If you really wanted to lean in to that aspect, you could sort of justify the zero-sum view, and maybe look for a way to sort and order by “how much” a team “won” each category. But again, what’s the right denominator for those measures? For instance, do we care about a team’s share of all Defensive Rebounds in the game? Or do we care about the share of Defensive Rebounds a team won relative to every opportunity it had to make a Defensive Rebound? How meaningful is ordering our rows by those kinds of shares? Even worse, some measures (notably Fouls) are bad to “win”, so we’d have to do something about those.

Our fundamental problem is that we just have two cases (the teams) and fifteen different measures, or variables. Each variable, except for the three percentages, is in effect on its own scale. There’s no direct way to make comparisons across them. Sure, some of these measures are probably going to be associated with one another—e.g. Turnovers and Points Off Turnovers—but the numeric values aren’t directly comparable in general. If you know a lot about basketball you might have some informative rules of thumb about each one of these measures, or some of them in combination. But at that point the lines in this particular graph are not going to be doing any work for you; you’ll just end up looking directly at the numbers. If we had data on all these measures for every NBA game for a whole season then we could of course do much more with them, because then each measure would have a distribution across all games and across all teams.

As it is, the purpose of the “Stats” screen in Apple Sports is just to summarize information from a single game. The other thing I could think of to do with the numbers as kind of graph is something like this:

This is marginally more helpful than the one before just because, again, it gets rid of the unhelpful zero-sum look of the original. As I hope you can immediately see, it creates many other difficulties. It also doesn’t do away with the core problem. That problem is principally one of information design rather than data visualization. What I mean is that what we’re trying to organize is, in effect, fifteen pairs of related but fundamentally distinct numbers. If we had fifteen cases and two variables things would be simple. But with fifteen variables and two cases … well, this is not the kind of thing you can make a single effective and non-confusing graph out of. That’s why I kind of sympathize with the designer. In a constrained space they have to show thirty numbers (thirty two, including the score). Lots of information. A straight table seems like it would be boring. Surely there’s some way to thematically integrate the numbers in a visually appealing manner that brings out some of the relationships across the rows. That’s what graphs do; it seems like the right thing to reach for. But at its heart this information is not a graph. It just sort of looks like one, and that ends up confusing people.

A few days ago, I presented Conformalized TabPFN: Prediction Intervals for a Pretrained Transformer for Tabular Data in Python and R. Today, it’s about TabICL, another state-of-the-art tabular foundation model. TabICL requires no token, as you’ll notice in the following Python and R code.

1 – Python version

!pip install tabicl nnetsauce # scikit-learn matplotlib numpy

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.linear_model import RidgeCV
from sklearn.metrics import mean_squared_error
from tabicl import TabICLRegressor
import nnetsauce as ns
import numpy as np
import matplotlib.pyplot as plt
from time import time

# ── data ───────────────────────────────────────────────────
X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# ── base models ────────────────────────────────────────────
models = {
    "TabICL": TabICLRegressor(),
    "RidgeCV": RidgeCV(),
}

results = {}
for name, reg in models.items():
    start = time()
    conf = ns.PredictionInterval(reg, level=95)
    conf.fit(X_train, y_train)
    pi = conf.predict(X_test, return_pi=True)
    print(f"{name:10s}  time={time() - start:.1f}s")

    coverage = np.mean((pi.lower <= y_test) & (pi.upper >= y_test))
    width    = np.mean(pi.upper - pi.lower)
    rmse     = np.sqrt(mean_squared_error(y_test, pi.mean))

    results[name] = {"pi": pi, "coverage": coverage,
                     "width": width, "rmse": rmse}
    print(f"{name:10s}  RMSE={rmse:.1f}  "
          f"coverage={coverage:.3f}  avg_width={width:.1f}")

# ── plot side-by-side ──────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
colors = {"TabICL": "orange", "RidgeCV": "steelblue"}
max_idx = 50

for ax, (name, res) in zip(axes, results.items()):
    pi = res["pi"]
    x  = range(max_idx)
    ax.fill_between(x, pi.lower[:max_idx], pi.upper[:max_idx],
                     alpha=0.35, color=colors[name], label="95% PI")
    ax.plot(x, pi.mean[:max_idx], "k--", lw=1.5, label="predicted")
    ax.plot(x, y_test[:max_idx], "k.", ms=6, alpha=0.4, label="observed")
    ax.set_title(
        f"{name}  |  cov={res['coverage']:.3f}  width={res['width']:.1f}"
    )
    ax.legend(fontsize=8)

plt.suptitle("Conformalized TabICL vs RidgeCV — diabetes dataset")
plt.tight_layout()
plt.show()

Checkpoint 'tabicl-regressor-v2-20260212.ckpt' not cached.
 Downloading from Hugging Face Hub (jingang/TabICL).




tabicl-regressor-v2-20260212.ckpt:   0%|          | 0.00/114M [00:00<?, ?B/s]


TabICL      time=21.8s
TabICL      RMSE=54.4  coverage=0.955  avg_width=226.1
RidgeCV     time=0.0s
RidgeCV     RMSE=53.9  coverage=0.955  avg_width=211.5

2 - R version

 %load_ext rpy2.ipython # in a Colab notebook, use this

%R install.packages("reticulate")

%%R  # in Colab/Jupyter with rpy2; remove this line for pure R

library(reticulate)

# pip install tabicl nnetsauce scikit-learn matplotlib numpy

sklearn_ds  <- import("sklearn.datasets")
sklearn_ms  <- import("sklearn.model_selection")
sklearn_m   <- import("sklearn.metrics")
sklearn_lm  <- import("sklearn.linear_model")
tabicl      <- import("tabicl")
ns          <- import("nnetsauce")
np          <- import("numpy")
plt         <- import("matplotlib.pyplot")

# ── data ───────────────────────────────────────────────────
d       <- sklearn_ds$load_diabetes(return_X_y = TRUE)
X <- d[[1]]; y <- d[[2]]
sp      <- sklearn_ms$train_test_split(X, y,
             test_size = 0.2, random_state = 42L)
X_train <- sp[[1]]; X_test <- sp[[2]]
y_train <- sp[[3]]; y_test <- sp[[4]]

# ── helper: fit + evaluate ─────────────────────────────────
eval_model <- function(reg, name) {
  conf <- ns$PredictionInterval(reg, level = 95L)
  conf$fit(X_train, y_train)
  pi   <- conf$predict(X_test, return_pi = TRUE)

  cov  <- np$mean((pi$lower <= y_test) * (pi$upper >= y_test))
  wid  <- np$mean(pi$upper - pi$lower)
  rmse <- sqrt(sklearn_m$mean_squared_error(y_test, pi$mean))

  cat(sprintf("%-10s  RMSE=%.1f  coverage=%.3f  avg_width=%.1f\n",
              name, rmse, cov, wid))
  invisible(pi)
}

# ── run both models ────────────────────────────────────────
pi_tabicl  <- eval_model(tabicl$TabICLRegressor(),  "TabICL")
pi_ridge   <- eval_model(sklearn_lm$RidgeCV(),       "RidgeCV")

# ── plot ───────────────────────────────────────────────────
max_idx <- 50L
x_range <- np$array(0:(max_idx - 1))

plot_pi <- function(pi, title, col) {
  x_fill <- np$concatenate(list(x_range, x_range[max_idx:1]))
  y_fill <- np$concatenate(list(
    pi$upper[1:max_idx], pi$lower[max_idx:1]))
  plt$fill(x_fill, y_fill, alpha=0.35, fc=col, ec="None", label="95% PI")
  plt$plot(x_range, pi$mean[1:max_idx], "k--", lw=1.5, label="predicted")
  plt$plot(x_range, y_test[1:max_idx], "k.", ms=6L, alpha=0.4, label="observed")
  plt$title(title); plt$legend(fontsize=8L)
}

fig <- plt$figure(figsize=c(12, 4))
plt$subplot(1L, 2L, 1L); plot_pi(pi_tabicl, "Conformalized TabICL", "orange")
plt$subplot(1L, 2L, 2L); plot_pi(pi_ridge,  "Conformalized RidgeCV", "steelblue")
plt$suptitle("Conformalized TabICL vs RidgeCV — diabetes dataset")
plt$tight_layout()
plt$show()

    WARNING: The R package "reticulate" only fixed recently
    an issue that caused a segfault when used with rpy2:
    https://github.com/rstudio/reticulate/pull/1188
    Make sure that you use a version of that package that includes
    the fix.
    TabICL      RMSE=54.4  coverage=0.955  avg_width=226.1
RidgeCV     RMSE=53.9  coverage=0.955  avg_width=211.5

Probably a dataset that’s too easy for a Transformer. Conformalizing simple models helps them, in general, to obtain coverage rates close to the nominal level, as we see for RidgeCV here.

The manifold hypothesis, the idea that real-world high-dimensional data concentrates near a low-dimensional curved subspace, is foundational to modern machine learning. Many popular manifold learning methods such as UMAP, t-SNE, Isomap, and diffusion maps do achieve dimensionality reduction by embedding data into a flat Euclidean space , but they do not attempt to directly learn the underlying manifold. In contrast, the 2025 paper by Robinett et al., Atlas-based manifold representations for interpretable Riemannian machine learning, offers a proof of concept for directly tackling the manifold hypothesis based on fundamental ideas from differential geometry. It provides an algorithm for learning a low dimensional manifold from point data by constructing an atlas of charts. The paper is also notable for the design of an efficient data structure for working with the learned atlas and for the extensive supplementary materials that include a GitHub Repository containing several practical Python algorithms for doing calculations on manifolds, and an extraordinary amount of implementation detail.

Reading through Robinett et al., however, requires a fairly deep background in the theory of differential geometry. This post is an attempt to provide an on-ramp to Robinett et al. by discussing the relatively simple example of the two dimensional sphere, embedded in . It implements the Atlas-Learn data structures and algorithms in R, uses them to learn and then goes on to validate the Atlas-Learn algorithm for the sphere via three independent methods: 1) use numerical integration along the manifold to trace a great circle on the sphere, 3) recover the radius of curvature of the sphere from the atlas, and 4) verify the Gauss-Bonnet Theorem for the sphere.

The R code was mostly worked out by Claude Sonnet 4.3 in the context of participating in the Posit beta test for its AI Assistant. I found the integration of the AI engine into the RStudio IDE an effective means of communicating with Claude and managing the project workflow.

library(tidyverse)
library(cluster)    # pam() for k-medoids partitioning
library(RANN)       # nn2() for fast k-nearest-neighbor queries
library(plotly)     # interactive 3D visualization
library(purrr)      # map() / imap() for list operations

#| label: atlas-functions

# ===========================================================================
# PART 1: Quadratic feature helpers
# ---------------------------------------------------------------------------
# These functions implement the d=2 specialization of the general quadratic
# feature map. For general d, phi(xi) would have choose(d+1, 2) components.
# For d=2 it has exactly 3: [xi1^2, xi1*xi2, xi2^2].
# ===========================================================================

# Maps xi in R^2 to the three quadratic monomials used to model surface curvature.
# General d would give choose(d+1,2) monomials; for d=2 this is exactly 3.
quad_features <- function(xi) {
  c(xi[1]^2, xi[1] * xi[2], xi[2]^2)
}

# Jacobian of quad_features w.r.t. xi: a 3x2 matrix (d=2 specialization).
# Row i is the gradient of the i-th monomial:
#   d/dxi [xi1^2]     = [2*xi1,   0   ]
#   d/dxi [xi1*xi2]   = [xi2,     xi1 ]
#   d/dxi [xi2^2]     = [0,       2*xi2]
quad_jacobian <- function(xi) {
  matrix(
    c(2 * xi[1],   xi[2],         0,
              0,   xi[1],  2 * xi[2]),
    nrow = 3, ncol = 2, byrow = TRUE
  )
}

# ===========================================================================
# PART 2: Core chart operations
# ---------------------------------------------------------------------------
# Each chart stores:
#   mean  : R^3 — chart center (centroid of the cluster)
#   L     : R^{3x2} — orthonormal tangent basis (columns = v1, v2 from SVD)
#   M     : R^{3x1} — unit surface normal (v3 from SVD; scalar normal because D-d=1)
#   K     : R^{1x3} — quadratic curvature coefficients (1 row because D-d=1)
#   ell_A : R^{2x2} — ellipsoidal domain matrix (2x2 because d=2)
# ===========================================================================

# Evaluate the inverse chart map: xi in R^2  ->  point in R^3.
# Formula: x = mean + L*xi + M * (K * phi(xi))
#   - mean + L*xi : linear move along the tangent plane
#   - M*(K*phi)   : quadratic normal correction encoding surface curvature
chart_eval <- function(chart, xi) {
  q <- quad_features(xi)          # 3-vector: [xi1^2, xi1*xi2, xi2^2]
  as.vector(
    chart$mean +
    chart$L %*% xi +              # 3x2 * 2x1 = 3x1 tangent contribution
    chart$M * as.numeric(chart$K %*% q)   # 3x1 * scalar normal correction
  )
}

# Jacobian of the inverse chart map at xi: a 3x2 matrix.
# J(xi) = L + M * (K * dQ(xi))
#   dQ(xi) is the 3x2 Jacobian of the monomial map (quad_jacobian).
#   K * dQ  is 1x3 * 3x2 = 1x2; then M * (K*dQ) is 3x1 * 1x2 = 3x2.
# This is the key object for geodesic integration: it maps tangent vectors
# in R^2 chart coordinates to ambient R^3 velocity vectors.
chart_jacobian <- function(chart, xi) {
  dQ <- quad_jacobian(xi)                   # 3x2 monomial Jacobian
  chart$L + chart$M %*% (chart$K %*% dQ)   # 3x2 total Jacobian J(xi)
}

# Project an ambient R^3 point onto this chart's tangent-plane coordinates (R^2).
# The linear projection xi = L^T * (x - mean) is exact for points on the tangent
# plane; for points on the actual surface it is a first-order approximation.
chart_project <- function(chart, x) {
  as.vector(t(chart$L) %*% (x - chart$mean))
}

# Test whether tangent coords xi lie within the chart's ellipsoidal domain.
# The domain is {xi in R^2 : xi^T * A * xi <= 1}, a 2D Mahalanobis ball.
# This is O(d^2) = O(4) for d=2.
chart_in_domain <- function(chart, xi) {
  as.numeric(t(xi) %*% chart$ell_A %*% xi) <= 1
}

# ===========================================================================
# PART 3: atlas_learn() — the main learning function
# ---------------------------------------------------------------------------
# Input : X (N x 3 matrix of surface points), k (number of charts)
# Output: an S3 object of class "atlas" containing k atlas_chart objects
#
# The algorithm runs four steps per chart:
#   (a) k-medoids partitioning
#   (b) local PCA to find the tangent plane and normal
#   (c) quadratic regression for the curvature coefficients K
#   (d) ellipsoidal domain construction
#
# d=2, D=3 specializations that are hard-coded here:
#   - SVD takes nv=3 (the full 3x3 right-singular-vector matrix)
#   - L = V[,1:2] is 3x2; M = V[,3] is 3x1 (one normal vector, not a matrix)
#   - nu (normal offset) is a scalar per point, not a vector
#   - K is fitted as a 1x3 row vector (one row per normal direction)
#   - ell_A is 2x2 (domain is a 2D ellipse)
# ===========================================================================

atlas_learn <- function(X, k, ellipsoid_scale = 1.1) {
  # X              : N x 3 matrix of surface points
  # k              : number of charts
  # ellipsoid_scale: inflate domains by this factor so adjacent charts overlap

  message("Fitting k-medoids (k = ", k, ") ...")
  # PAM (Partitioning Around Medoids) is preferred over k-means for surfaces:
  # medoids are actual data points, making the partition robust to outliers.
  km <- cluster::pam(X, k, diss = FALSE)

  charts <- seq_len(k) |>
    purrr::map(\(j) {

      # --- (a) Extract the cluster and center it ----------------------------
      idx <- which(km$clustering == j)
      Xj  <- X[idx, , drop = FALSE]   # N_j x 3 matrix of cluster points
      m   <- colMeans(Xj)             # chart center (3-vector)
      Xc  <- sweep(Xj, 2, m)         # centered cluster: N_j x 3

      # --- (b) Local PCA: tangent plane and normal -------------------------
      # SVD of the centered cluster reveals local geometry:
      #   - first two right singular vectors (v1, v2) span the tangent plane
      #   - third right singular vector (v3) is the surface normal
      # d=2, D=3: we always take nv=3 because D=3 (fully determined).
      sv  <- svd(Xc, nu = 0, nv = 3)

      # L: 3x2 tangent basis (orthonormal columns). General form: R^{D x d}.
      L   <- sv$v[, 1:2, drop = FALSE]

      # M: 3x1 unit normal. d=2, D=3 specialization: D-d=1, so there is
      # exactly ONE normal direction and M is a column vector, not a matrix.
      M   <- sv$v[, 3, drop = FALSE]

      # Project each centered point into tangent / normal coordinates.
      # tau: N_j x 2 — tangent coordinates (2D because d=2)
      # nu:  N_j x 1 — normal offsets (scalar because D-d=1)
      tau <- Xc %*% L    # N_j x 2
      nu  <- Xc %*% M    # N_j x 1 (scalar per point; would be N_j x (D-d) in general)

      # --- (c) Quadratic curvature regression ------------------------------
      # Fit: nu ~ K * phi(tau)  where phi(tau) = [tau1^2, tau1*tau2, tau2^2]
      # K is 1x3 here (one output because D-d=1; would be (D-d)x3 in general).
      #
      # Design matrix Q: N_j x 3, each row is phi(tau_i)
      Q   <- t(apply(tau, 1, quad_features))   # N_j x 3

      # Ridge-regularized least squares: K = nu^T * Q * (Q^T*Q + esp*I)^{-1}
      # The ridge penalty eps=1e-10 prevents singularity when clusters are
      # nearly collinear in tangent coordinates.
      K   <- t(nu) %*% Q %*% solve(crossprod(Q) + 1e-10 * diag(3))  # 1 x 3

      # --- (d) Ellipsoidal domain ------------------------------------------
      # Domain: {xi in R^2 : xi^T A xi <= 1}, a Mahalanobis ball.
      # A_raw = Cov(tau)^{-1}: inverse covariance of the tangent coordinates.
      # This adapts the domain shape to the data spread in each direction.
      A_raw <- solve(cov(tau))

      # Scale A so the outermost point lands on the boundary, then inflate
      # by ellipsoid_scale (default 1.1) to ensure overlap with neighbors.
      qvals <- apply(tau, 1, \(xi) as.numeric(t(xi) %*% A_raw %*% xi))
      ell_A <- A_raw / (ellipsoid_scale * max(qvals))   # 2x2 positive-definite

      # Pack all chart data into an S3 object
      structure(
        list(mean = m, L = L, M = M, K = K, ell_A = ell_A, idx = idx,
             n_points = length(idx)),
        class = "atlas_chart"
      )
    })

  # Return the full atlas as an S3 object
  structure(
    list(
      charts        = charts,
      k             = k,
      clustering    = km$clustering,
      n_points      = nrow(X),
      ambient_dim   = ncol(X),   # D = 3
      intrinsic_dim = 2L         # d = 2 (hard-coded for surfaces in R^3)
    ),
    class = "atlas"
  )
}

# ===========================================================================
# PART 4: S3 print / summary methods
# ===========================================================================

print.atlas_chart <- function(x, ...) {
  cat(sprintf(
    "<atlas_chart>  %d points | mean [%s] | cond(ell_A) = %.1f\n",
    x$n_points,
    paste(round(x$mean, 3), collapse = ", "),
    kappa(x$ell_A)
  ))
  invisible(x)
}

# Returns a per-chart summary tibble
summary.atlas <- function(object, ...) {
  purrr::imap_dfr(object$charts, \(ch, i)
    tibble::tibble(
      chart      = i,
      n_points   = ch$n_points,
      mean_norm  = round(sqrt(sum(ch$mean^2)), 4),
      cond_ell_A = round(kappa(ch$ell_A), 1)
    )
  )
}

print.atlas <- function(x, ...) {
  cat(sprintf(
    "<atlas>  k = %d | ambient R^%d | intrinsic R^%d | %d points\n\n",
    x$k, x$ambient_dim, x$intrinsic_dim, x$n_points
  ))
  print(summary(x))
  invisible(x)
}

# ===========================================================================
# PART 5: Chart lookup
# ===========================================================================

# Return the index of the nearest chart whose domain contains x.
# Falls back to the globally nearest chart center if none qualify.
# Checking only the 6 nearest charts (by Euclidean center distance) is
# sufficient in practice and avoids an O(k) domain test at every step.
find_chart <- function(atlas, x) {
  dists      <- map_dbl(atlas$charts, \(ch) sum((x - ch$mean)^2))
  candidates <- order(dists)[seq_len(min(6L, atlas$k))]
  in_domain  <- purrr::keep(
    candidates,
    \(i) chart_in_domain(atlas$charts[[i]],
                         chart_project(atlas$charts[[i]], x))
  )
  if (length(in_domain) > 0L) in_domain[[1L]] else which.min(dists)
}

# ===========================================================================
# PART 6: Quasi-Euclidean retraction
# ---------------------------------------------------------------------------
# Advances one step from (chart_idx, xi) in the ambient direction tau_r3 (R^3).
#
# The pseudoinverse J^+ = (J^T J)^{-1} J^T is the key d=2, D=3 quantity:
#   - J is 3x2 (D x d)
#   - J^T J is 2x2 (d x d) — the pullback metric g; invertible for full-rank J
#   - J^+ is 2x3 (d x D) — pulls ambient vectors back to chart coordinates
#
# In general D > d, J^+ is the minimum-norm left inverse of J.
# ===========================================================================

retract_step <- function(atlas, chart_idx, xi, tau_r3) {
  chart  <- atlas$charts[[chart_idx]]
  J      <- chart_jacobian(chart, xi)        # 3x2 (D x d for d=2, D=3)
  Jp     <- solve(crossprod(J)) %*% t(J)    # 2x3 pseudoinverse: (J^T J)^{-1} J^T
  xi_new <- xi + as.vector(Jp %*% tau_r3)   # advance in chart coordinates

  if (chart_in_domain(chart, xi_new)) {
    # Still inside the same chart: evaluate and return
    return(list(
      chart_idx = chart_idx,
      xi        = xi_new,
      x         = chart_eval(chart, xi_new)
    ))
  }

  # Step crossed a chart boundary: find a new host chart and re-project
  x_cand <- chart_eval(chart, xi_new)       # approximate ambient position
  ci_new <- find_chart(atlas, x_cand)       # nearest chart that contains x_cand
  xi2    <- chart_project(atlas$charts[[ci_new]], x_cand)   # re-project
  list(
    chart_idx = ci_new,
    xi        = xi2,
    x         = chart_eval(atlas$charts[[ci_new]], xi2)
  )
}

# ===========================================================================
# PART 7: Geodesic path integration
# ---------------------------------------------------------------------------
# Traces a geodesic from x_start in direction direction_r3 (ambient R^3).
# Returns a tibble with columns x, y, z, step, chart_idx.
#
# At each step the ambient velocity vector is re-projected into the current
# chart's tangent plane and pushed back to ambient.  This keeps the direction
# of motion consistent across chart boundaries — without it, the fixed tau
# vector drifts after every chart transition and the path wanders off the
# intended geodesic.
# ===========================================================================

atlas_geodesic <- function(atlas, x_start, direction_r3,
                           n_steps = 100L, step_size = 0.02) {
  ci  <- find_chart(atlas, x_start)
  ch  <- atlas$charts[[ci]]
  xi  <- chart_project(ch, x_start)

  # Project the initial direction onto the starting chart's tangent plane,
  # then normalize to unit speed.  This ensures tau_r3 is always in T_x M.
  J      <- chart_jacobian(ch, xi)
  Jp     <- solve(crossprod(J)) %*% t(J)          # 2x3 pseudoinverse
  tau_r3 <- as.vector(J %*% (Jp %*% direction_r3))  # project onto tangent plane
  tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size  # unit-speed scaling

  steps     <- vector("list", n_steps + 1L)
  chart_ids <- integer(n_steps + 1L)
  steps[[1L]]     <- chart_eval(ch, xi)
  chart_ids[[1L]] <- ci

  for (i in seq_len(n_steps)) {
    res    <- retract_step(atlas, ci, xi, tau_r3)
    ci     <- res$chart_idx
    xi     <- res$xi
    steps[[i + 1L]]     <- res$x
    chart_ids[[i + 1L]] <- ci

    # Re-project tau_r3 into the current chart's tangent plane and normalize.
    # This is the identity-transport approximation: it holds the direction
    # constant in the current chart's frame rather than applying the
    # Christoffel correction that true parallel transport would require.
    J_new  <- chart_jacobian(atlas$charts[[ci]], xi)
    Jp_new <- solve(crossprod(J_new)) %*% t(J_new)
    tau_r3 <- as.vector(J_new %*% (Jp_new %*% tau_r3))
    tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size
    J      <- J_new
  }

  do.call(rbind, steps) |>
    as_tibble(.name_repair = "minimal") |>
    set_names(c("x", "y", "z")) |>
    mutate(step = row_number(), chart_idx = chart_ids)
}

# ===========================================================================
# PART 8: Sphere-specific helpers
# ---------------------------------------------------------------------------
# These are ground-truth functions for S^2 used only for validation;
# they play no role in the Atlas-Learn algorithm itself.
# ===========================================================================

# Great-circle geodesic distance between two unit vectors (in radians).
# d(p, q) = arccos(p . q), clamped to [-1, 1] to avoid NaN from floating point.
sphere_dist <- function(p, q) {
  acos(pmax(pmin(sum(p * q), 1.0), -1.0))
}

# Unit tangent vector at p pointing toward q along the great circle.
# Computed by projecting q onto the tangent plane at p and normalizing.
geodesic_direction_sphere <- function(p, q) {
  v   <- q - sum(p * q) * p   # remove radial component
  nrm <- sqrt(sum(v^2))
  if (nrm < 1e-12) return(NULL)
  v / nrm
}

# Run the Atlas geodesic from p toward q for the exact number of steps
# needed to cover the true great-circle distance, then measure endpoint error.
atlas_endpoint_error <- function(atlas, p, q, step_size = 0.02) {
  d   <- sphere_dist(p, q)
  dir <- geodesic_direction_sphere(p, q)

  if (d < 1e-6 || is.null(dir)) {
    return(tibble(true_dist = d, endpoint_error = 0, n_steps = 0L))
  }

  n_steps  <- max(1L, round(d / step_size))
  path     <- atlas_geodesic(atlas, p, dir,
                              n_steps = n_steps, step_size = step_size)

  endpoint <- as.numeric(path[nrow(path), c("x", "y", "z")])
  # Re-project endpoint onto S^2 before measuring angle error.
  # This separates geodesic-direction error from on-manifold drift.
  endpoint <- endpoint / sqrt(sum(endpoint^2))

  tibble(
    true_dist      = d,
    endpoint_error = sphere_dist(endpoint, q),
    n_steps        = as.integer(n_steps)
  )
}

#| code-summary: "Show the code"

#| label: generate-sphere

set.seed(4471)
N     <- 2000L
theta <- runif(N, 0, 2 * pi)
phi   <- acos(runif(N, -1, 1))   # uniform area measure on S^2

sphere_pts <- tibble(
  x = sin(phi) * cos(theta),
  y = sin(phi) * sin(theta),
  z = cos(phi)
)

X_sphere <- as.matrix(sphere_pts)

#| label: learn-atlas

atlas_cache <- "atlas_s2_cache.rds"

if (file.exists(atlas_cache)) {
  message("Loading cached atlas from ", atlas_cache)
  atlas_s2 <- readRDS(atlas_cache)
} else {
  message("Cache not found — fitting atlas (k = 25) ...")
  atlas_s2 <- atlas_learn(X_sphere, k = 25L)
  saveRDS(atlas_s2, atlas_cache)
  message("Atlas saved to ", atlas_cache)
}

# Attach chart labels to the point-cloud tibble for visualization
sphere_pts <- sphere_pts |>
  mutate(chart = factor(atlas_s2$clustering))

atlas_s2
├── k              integer          number of charts (25)
├── n_points       integer          total data points (2000)
├── ambient_dim    integer          dimension of embedding space (3 = D)
├── intrinsic_dim  integer          dimension of the manifold (2 = d)
├── clustering     integer[2000]    chart label for every data point
└── charts         list[25]         one atlas_chart per cluster
    └── charts[[j]]
        ├── mean     numeric[3]     centroid of cluster j in R³
        ├── L        numeric[3×2]   orthonormal tangent basis (columns = v₁, v₂)
        ├── M        numeric[3×1]   unit surface normal = v₃  [For d=2: single vector]
        ├── K        numeric[1×3]   quadratic curvature coefficients  [for d=2: 1 row]
        ├── ell_A    numeric[2×2]   ellipsoidal domain matrix  [d=2 specialisation: 2×2]
        ├── idx      integer[n_j]   row indices into the N×3 data matrix
        └── n_points integer        cluster size n_j

purrr::imap_dfr(atlas_s2$charts, function(ch, j) {
  tibble(
    Chart  = j,
    `n pts` = ch$n_points,
    `centre (x,y,z)` = paste(round(ch$mean, 2), collapse = ", "),
    # For the unit sphere, M should align with the radial direction (M·r̂ ≈ 1)
    `M · r̂`  = round(abs(sum(ch$M * ch$mean / sqrt(sum(ch$mean^2)))), 3),
    # Quadratic curvature coefficients
    k1 = round(ch$K[1], 3),
    k2 = round(ch$K[2], 3),
    k3 = round(ch$K[3], 3),
    # Semi-axes of the 2D ellipsoidal domain (1/sqrt of eigenvalues of ell_A)
    `ell axis a` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[2]), 3),
    `ell axis b` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[1]), 3)
  )
}) |> print(n = 25)

# A tibble: 25 × 9
   Chart `n pts` `centre (x,y,z)`    `M · r̂`     k1     k2     k3 `ell axis a`
   <int>   <int> <chr>                 <dbl>  <dbl>  <dbl>  <dbl>        <dbl>
 1     1      82 0.91, 0.3, -0.16      0.999 -0.216 -0.078 -0.068        0.498
 2     2      70 0.35, -0.9, -0.06     1     -0.136  0     -0.163        0.486
 3     3      84 -0.64, 0.71, -0.13    1      0.182 -0.031  0.101        0.44 
 4     4      73 -0.04, 0.95, -0.12    1      0.162 -0.024  0.112        0.463
 5     5      81 -0.75, 0.4, 0.46      1      0.115 -0.056  0.184        0.431
 6     6      73 0.28, 0.6, -0.7       0.999 -0.203  0.016 -0.114        0.518
 7     7      78 -0.93, 0.07, -0.22    1      0.169 -0.087  0.055        0.46 
 8     8      86 -0.08, -0.3, 0.91     1     -0.237  0.051  0.099        0.53 
 9     9      73 -0.18, -0.87, -0.38   1      0.111 -0.036  0.127        0.421
10    10      85 0.87, 0.07, 0.4       1     -0.148 -0.002 -0.13         0.471
11    11      75 -0.36, 0.21, 0.87     1      0.2   -0.034  0.054        0.556
12    12      67 -0.79, -0.52, -0.21   1      0.178 -0.001  0.109        0.431
13    13      84 0.04, 0.05, -0.96     1      0.182  0.012  0.081        0.446
14    14      66 -0.39, -0.83, 0.31    1     -0.189 -0.006 -0.066        0.511
15    15      84 0.28, -0.55, -0.74    1      0.132 -0.066  0.105        0.431
16    16      87 -0.52, 0.43, -0.68    1     -0.251  0.009  0.026        0.511
17    17     100 0.83, -0.44, -0.2     1      0.144  0.082  0.133        0.526
18    18      68 0.07, -0.8, 0.55      1     -0.24   0.072 -0.088        0.474
19    19      91 0.59, 0.71, 0.24      0.999  0.229 -0.002 -0.003        0.537
20    20      93 -0.81, -0.26, 0.45    1      0.17   0.084  0.115        0.455
21    21      95 0.41, 0.27, 0.83      1      0.155 -0.087  0.056        0.532
22    22      66 -0.46, -0.32, -0.78   1     -0.179 -0.025 -0.072        0.525
23    23      70 -0.12, 0.78, 0.54     1      0.183  0.045  0.1          0.495
24    24      93 0.67, 0.07, -0.69     1      0.192 -0.001  0.116        0.453
25    25      76 0.62, -0.45, 0.58     1     -0.215 -0.039  0.028        0.512
# ℹ 1 more variable: `ell axis b` <dbl>

#| label: atlas-chart-plot

n_charts <- length(unique(sphere_pts$chart))
chart_colors <- colorRampPalette(RColorBrewer::brewer.pal(9, "Set1"))(n_charts)

plot_ly(
  data   = sphere_pts,
  x = ~x, y = ~y, z = ~z,
  type   = "scatter3d",
  mode   = "markers",
  color  = ~chart,
  colors = chart_colors,
  marker = list(size = 2.5, opacity = 0.6),
  text   = ~paste0("Chart ", chart),
  hoverinfo = "text"
) |>
  layout(
    title  = "S² coloured by Atlas-Learn chart membership (k = 25)",
    scene  = list(
      xaxis = list(title = "x"),
      yaxis = list(title = "y"),
      zaxis = list(title = "z"),
      aspectmode = "cube"
    ),
    showlegend = FALSE
  )

Here I am again, writing about crochet and programming! I’ve continued creating my tons of shitty stitches cute creatures.

New Git analogy! The crochet lifeline

First Git/crochet analogy.

I read a great crochet book about creating your own amigurumi patterns: Créer ses propres modèles d’amigurumis au crochet by Clotilde Massot and Lise Grandjonc. One of its authors (Clotilde Massot) is a software developer who, among other things, published an octocat pattern that I bought and used… Anyway in the book they explain that when you create a pattern, you will probably have to undo your work several times. Undoing a round is easy if you have a stitch marker in the first stitch of the current/last round: you undo until you hit that stitch marker. But what about undoing several rounds? In that case, you’ll be better off if instead of using a stitch marker at the beginning of the current round, you use contrasting yarn stuck under the first stitch of each round. That yarn creates what the authors of the book call your lifeline!

Now, if that’s not a great analogy for commits and the ability to reset…

Communities of practice

If you’ve had the opportunity to attend or watch the wonderful useR! 2025 keynote talk of my rOpenSci colleague Yanina Bellini Saibene, you’ve heard of communities of practice. In her talk, Yani quoted Etienne Wenger who defined communities of practice as “groups of people who share a passion for something that they know how to do, and who interact regularly in order to learn how to do it better“. Yani mentioned her swimming team and English conversation club. Well, I found a community of practice for crochet: a stitch club at a local café! Participants meet up to crochet side by side, talking a lot about crochet: comparing yarns, exchanging tips, etc. The first time I went, we even started with a round of introductions where we said what each of us would work on during that meeting, which reminded me of rOpenSci coworking sessions where participants do exactly that.

Usefulness of seeing others’ work

I’ve worked on some patterns by Yan Schenkel a.k.a. Pica Pau. One cool aspect of the patterns is that they include a link to a gallery where anyone can upload pics of their take on each creature. So you can look at them, maybe seeing a crucial (to you) angle that’s absent from the pattern, noticing whether some “flaw” of your own project is present in others’ projects, comparing variations and picking what you prefer before you start, etc. For instance I stared at many pics of Alberto Seagull before making the legs for mine.

The usefulness of seeing others’ crocheted animals reminds me of how reading and reviewing open-source, or our colleagues’, code helps us (and LLMs, I suppose) learn how to do, or not to do, some things, how it help us refine our taste! Ironically, I do publish most of my code, but I haven’t had interest in doing that for my crocheting yet.

Conclusion

In summary, I keep finding excuses to talk about crochet. My pink octocat approves!

Note for myself and others: how to remove the “Leaflet | ” prefix in map attribution using the R package {leaflet}.

According to some sources we could write leaflet(options = leafletOptions(attributionPrefix = "")) but it doesn’t work in my case.

So instead we can execute some javascript:

library(leaflet)

leaflet() |>
  addTiles(urlTemplate = "", attribution = "Only my data") |> 
  htmlwidgets::onRender("function(el, x) {
    // Remove Leaflet attribution prefix
    this.attributionControl.setPrefix('');
    }")

Introduction

rOpenSci’s R-universe system is an open source platform allowing users to create their own CRAN-like universe of R packages.

It is absolutely fantastic. It is particularly useful in one area I research, Mendelian randomization (at the interface of Epidemiology and Genetic Epidemiology), because a lot of the packages are GitHub/GitLab-only.

Therefore, I setup and maintain https://mrcieu.r-universe.dev/ to include both packages from our MRCIEU GitHub organisation (from the MRC Integrative Epidemiology Unit at the University of Bristol, UK), and as many of the GitHub-only packages for Mendelian randomization I could find.

It is difficult to overstate how useful this is. For the first time, not only do researchers have a list of the Mendelian randomization packages in one place but they can install binaries – without having to go through the hassle – especially on (Ubuntu) Linux – of remotes::install_github(). Researchers can also see how often packages are updated and R-universe checks for changes in packages approximately every hour, keeping it always up to date.

This post gives five tips I have developed to help manage my R-universe.

Tip 1: Referring to a package from a pull request instead of from a branch on a fork

In the Mendelian randomization field many of these GitHub-only packages are not well written or abandoned once the PhD student/researcher leaves. Often when I add a package to our R-universe I find that their build fails, or they have R CMD check errors and warnings, or after several months their build fails because they are not maintained. I sometimes look into the failed builds and check problems. If it’s clear just a few fixes are required to rectify the situation I often open a pull request. Often that pull request is not responded to.

Previously, for such cases I would switch the source of the package entry in packages.json to be from the relevant branch on my fork. However, I have always felt a bit uneasy about this. I wondered if GitHub had a way to refer to the pull request branch without having to switch the repository. It turns out that it does. The format of pull request branch names is refs/pull/{number}/head where {number} is the number assigned once the PR is opened. Therefore, when I open a PR on a package I now add the "branch" field to the package entry in packages.json as follows.

  {
    "package": "GWASBrewer",
    "url": "https://github.com/jean997/GWASBrewer",
    "branch": "refs/pull/18/head"
  },

I switch back to the default branch if the PR is merged.

Tip 2: Justfile recipe for adding a package to packages.json

I regularly find that I need to add or remove a package. Manually editing the packages.json file is not hard, but I have found the following Justfile (Just is like Make, but specifically designed for running commands and has a much friendlier syntax) recipes helpful for doing this quickly.

These recipes require uv and just to be installed and on your PATH (uv automatically installs the required version of Python and creates/destroys/manages any required virtual environments). To use them, copy them into a text file named justfile at the top level of your R-universe registry repository and follow the instructions.

This recipe adds a package to your packages.json in alphabetical order. It has one required argument and 3 optional arguments.

# add a package entry to packages.json in alphabetical order
[arg("branch", short="b")]
[arg("pkgname", short="p")]
[arg("subdir", short="s")]
add url pkgname="" branch="" subdir="":
    #!/usr/bin/env -S uv run --python 3.14 python3
    import json, re, sys
    url = "{{ url }}"
    if re.fullmatch(r'[^/]+/[^/]+', url):
        url = f"https://github.com/{url}"
    pkgname = "{{ pkgname }}" or url.rstrip("/").split("/")[-1]
    branch = "{{ branch }}"
    subdir = "{{ subdir }}"
    with open("packages.json") as f:
        packages = json.load(f)
    if any(p["package"] == pkgname for p in packages):
        print(f"Error: '{pkgname}' already exists in packages.json", file=sys.stderr)
        sys.exit(1)
    entry = {"package": pkgname, "url": url}
    if branch:
        entry["branch"] = branch
    if subdir:
        entry["subdir"] = subdir
    packages.append(entry)
    packages.sort(key=lambda p: p["package"].lower())
    with open("packages.json", "w") as f:
        json.dump(packages, f, indent=2)
        f.write("\n")
    print(f"Added {pkgname}")

Where url is say https://github.com/MRCIEU/TwoSampleMR, except that for GitHub packages you can specify this as MRCIEU/TwoSampleMR.

To add a GitHub package whose name matches its repository name, simply run

just add username/reponame

You can inspect the recipe’s arguments and options with

just --usage add

Usage: just add [OPTIONS] url

Arguments:
  url

Options:
  -p pkgname [default: ""]
  -b branch [default: ""]
  -s subdir [default: ""]

The 3 optional arguments allow you to specify the package name (-p pkgname), branch (-b branchname), or subdirectory (-s subdirectory) the package is in. For example, to add a GitHub package whose package name does not match its repository name run

just add username/reponame -p pkgname

Tip 3: Justfile recipe for removing a package from packages.json

This recipe removes a package from your packages.json.

# remove a package entry from packages.json
remove pkgname:
    #!/usr/bin/env -S uv run --python 3.14 python3
    import json, sys
    pkgname = "{{ pkgname }}"
    with open("packages.json") as f:
        packages = json.load(f)
    filtered = [p for p in packages if p["package"] != pkgname]
    if len(filtered) == len(packages):
        print(f"Error: '{pkgname}' not found in packages.json", file=sys.stderr)
        sys.exit(1)
    with open("packages.json", "w") as f:
        json.dump(filtered, f, indent=2)
        f.write("\n")
    print(f"Removed {pkgname}")

Run it with

just remove pkgname

Tip 4: Justfile recipe for checking packages.json is valid

When manually editing packages.json it is very easy to forget a comma or to miss a closing bracket or quotation mark. This recipe checks your JSON is valid.

# check packages.json
check:
    uv run --python 3.14 -m json.tool packages.json > /dev/null && echo "JSON check passed"

Run it with

just check

Tip 5: Conveniently view a package’s dependencies

Knowing a package’s full strong dependency list is useful — for example, when a breaking change somewhere in the chain causes unexpected build failures. While there are several ways to determine this in R, R-universe shows you the full list immediately.

Navigate to the R-universe page for the package you are interested in, say https://mrcieu.r-universe.dev/TwoSampleMR and click the dependencies pill.

It expands showing the full dependency list.

Summary

In summary, I have shown five tips I find useful to manage a large R-universe.